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Abstract 

We study the dynamic critical behavior of a Swendsen-Wang-type algo- 
rithm for the Ashkin-Teller model. We find that the Li-Sokal bound on the 
autocorrelation time (r m t t £ > const X C'h) holds along the self-dual curve of 
the symmetric Ashkin-Teller model, and is almost but not quite sharp. The 
ratio 7i ntj ,f /(7f/ appears to tend to infinity either as a logarithm or as a small 
power (0.05 < p < 0.12). In an appendix we discuss the problem of extracting 
estimates of the exponential autocorrelation time. 

Key Words: Ashkin-Teller model, Ising model, Potts model, Monte Carlo, 
dynamical critical behavior, cluster algorithm, Swendsen-Wang algorithm, Li- 
Sokal bound, critical slowing-down, autocorrelation time, fitting correlated 
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1 Introduction 



Critical slowing-down has become one of the main limiting factors of the state of 
art of Monte Carlo simulations [1,2,3,4]. The autocorrelation time r, which roughly 
measures the Monte Carlo time between two statistically independent configura- 
tions, diverges near a critical point. More precisely, for a finite system of linear 
size L at criticality, we expect a behavior r ~ L z for large L: here the power z 
is a dynamic critical exponent, which characterizes the dynamic universality class 
of the Monte Carlo algorithm. The traditional local algorithms (such as single- 
site Metropolis) have a dynamic critical exponent z > 2. This is a severe critical 
slowing-down, in which the amount of computer work needed to study a lattice of 
size L grows approximately L 2 times faster than the naive geometrical factor L d 
(d being the dimensionality of the lattice). To study the static critical behavior 
we need high-precision data (run lengths > 10 3 r). In practice, it is very difficult 
to obtain high-precision data for large lattices with this kind of algorithm. To 
study the dynamic critical behavior, the situation is even worse, as we need much 
higher statistics (run lengths J> 10 4 r). The geometrical factor L d is unavoidable for 
the usual Monte Carlo simulations, so the elimination (or reduction) of the critical 
slowing-down is the only way to make Monte Carlo simulations feasible close to a 
critical point. 

Dramatic progress in this direction was stimulated by the introduction of the 
so-called cluster algorithms [5]. Instead of sequentially updating the whole lattice 
by single-spin moves, these algorithms employ non-local moves, such as cluster 
flips. For the ferromagnetic g-state Potts model, the Swendsen-Wang (SW) cluster 
algorithm [5] achieves a a significant reduction in z compared to the local algorithms: 
one has z between and 1, where the exact value depends on the number of Potts 
states and on the dimensionality of the lattice [4]. The two-dimensional (2D) Ising 
model is the most favorable case: the critical slowing-down becomes extremely 
weak, with estimates from different workers ranging from z = 0.35 ± 0.01 [5] to 
z = 0.25 ± 0.01 [6,7] to z = X log (i.e., r ~ logL) [8]. Unfortunately, it is very 
hard to distinguish between the power-law and logarithmic scenarios using only 
lattices with L < 512 [6,7]. In other cases, the performance of the SW algorithm 
is less impressive (though still quite good): e.g., z = 0.55 ± 0.03 for the 2D 3-state 
Potts model [9], and z ~ 1 for the 2D 4-state Potts model [9] and for the 4D Ising 
model [10,11]. Clearly, we would like to understand why this algorithm works so 
well in some cases and not in others; we hope in this way to obtain new insights 
into the dynamics of non-local Monte Carlo algorithms, with the ultimate aim of 
devising new and more efficient algorithms. 

A single-cluster variant of the SW algorithm was introduced by Wolff [12]. In- 
stead of updating all the clusters (with a given probability), only one cluster is 
selected and updated. It is not known why the dynamic exponents z\c associated 
to the single-cluster dynamics are very close to those of the SW dynamics in some 
cases (e.g., 2D q = 2,3 Potts models) but not in other cases (e.g., Ising model in 
dimension d > 3) [12,13]. A priori one would expect the two algorithms to belong 
to different dynamic universality classes. 
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There is at present no adequate theory for predicting the dynamic critical be- 
havior of an SW-type algorithm. However, there is one rigorous lower bound on 
z. In 1989 Li and Sokal [9] showed that the autocorrelation times of the standard 
(multi-cluster) SW algorithm for the ferromagnetic g-state Potts model are bounded 
below by a multiple of the specific heat: 

Tint,AA, T int:£ , r exp > const X C H • (1-1) 

Here Af is the bond density in the SW algorithm, £ is the energy, and Ch is the 
specific heat; T int and r exp denote the integrated and exponential autocorrelation 
times, respectively [1,4]. As a result one has 

a 

z iat,Afj z iat,£j z exp > ~ , (1-2) 

where a and v are the standard static critical exponents. Thus, the SW algorithm 
for the g-state Potts model cannot completely eliminate the critical slowing-down in 
any model in which the specific heat is divergent at criticality (although one might 
hope to obtain z = X log in the the 2D Ising model, where the specific heat is 
only logarithmically divergent). Now, one would like to know whether the bound 
(1.2) on critical exponents is sharp: that is, does it hold as equality, or only as a 
strict inequality? In more detail, one would like to know whether the bound (1.1) 
on the autocorrelation times is sharp (t/Ch bounded), sharp modulo a logarithm 
(r I Ch ~ log p L), or not sharp (r /Ch ~ L p with p > 0). 

Unfortunately, the empirical situation for the 2D Potts models is not very clear. 
For the Ising case, the bound (1.2) would be sharp if (and only if) the autocorrelation 
time grows like a logarithm; this is consistent with the data but not demanded by 
it [5,6,7,8]. For the 3-state Potts model, the bound is apparently not sharp: we 
have z = 0.55 ± 0.03 [9] versus a/v = | = 0.4 [14]. The 4-state Potts model is 
rather peculiar: the naive fit to the data, z = 0.89 ± 0.05 [9], is smaller than the 
(exactly known) value of a/v = 1 [15]. The explanation of this paradox is that the 
true leading term in the specific heat has a multiplicative logarithmic correction, 
C H ~ Llog- 3/2 L [16,17,18]; and indeed the observed exponent a/v (from a naive 
power-law fit) is 0.75 ± 0.01 [9], consistent with the bound (1.2). It is reasonable to 
conjecture that the true behavior of the autocorrelation time is likewise of the form 
t ~ Llog p L (with p > — |), in which case the bound (1.1) would be sharp modulo 
a multiplicative logarithm. 

So we are in a strange situation: the Li-Sokal bound might be sharp (possibly 
modulo a logarithm) for the 2D Potts models with q = 2 and q = 4, but it is 
apparently not sharp for q = 3. 1 

There is yet another way of "interpolating" between the 2-state (Ising) and 4- 
state Potts models: both are special cases of the Ashkin-Teller (AT) model [23], 

1 For Ising models in lattice dimensions d > 3, the bound (1.2) is clearly not sharp. For the 3D 
Ising model, estimates of z range from 0.339±0.004 to 0.75±0.01 [5,19,20], while a/v « 0.17 [21,22]. 
For Ising models in dimension d > 4, we expect z = 1 (possibly modulo a logarithm in d = 4) [10,11], 
while of course a jv = (or x log 2/3 in d = 4). 
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which has two interacting Ising spins on each lattice site. (For a review of the AT 
model, see Section 2 below.) The symmetric AT model (which is equivalent to the 
general Z4 clock model) presents a very rich phase diagram. In particular, one of the 
critical curves (namely, the self-dual curve) is quite unusual: the critical exponents 
vary continuously along this curve (as a result of a marginal operator), thus violating 
the usual notion of universality. One point on this critical curve is precisely the 4- 
state Potts model at criticality, while another point on this curve corresponds to a 
pair of decoupled Ising models. Thus, new insights on the sharpness of the Li-Sokal 
bound in the 2-state and 4-state Potts models might be obtainable by studying the 
same question on the self-dual curve of the symmetric AT model. 

Wiseman and Domany [24] devised the first SW-type algorithm for the AT 
model. Though their method of derivation is rather complicated, the algorithm is 
simple and reduces to the well-known SW algorithms in the special cases of the Ising 
and 4-state Potts models. The same algorithm had been independently introduced 
by Laanait et al. in another context [25]: they studied a model closely related to the 
AT model, and they used the same SW-type algorithm as a tool for their rigorous 
proofs. 

In this paper we would like to address the issue of the sharpness of the Li-Sokal 
bound (1.1)/ (1.2) along the self-dual curve of the symmetric AT model, and in 
particular for the 4-state Potts model. We propose two variants of the algorithm. 
The first, which we call the "direct" algorithm, is essentially the same as that of 
Wiseman and Domany [24]; however, we think that our derivation is simpler. (The 
reader can judge!) The second variant, which we call the "embedding" algorithm, 
is somewhat simpler to implement in practice; it is not equivalent to the direct 
algorithm, although we expect it to lie in the same dynamic universality class. 

We have studied numerically the multi-cluster ("standard SW") version of the 
embedding algorithm 2 at three points on the AT self-dual curve: the 4-state Potts 
model and two additional models (ZF and X2) interpolating between the 4-state 
Potts and Ising models. We have used lattices up to L = 512 (as well as L = 1024 
for the 4-state Potts model), and have systematically employed finite-size-scaling 
techniques to analyze the numerical data. We have also reanalyzed the very precise 
data reported by Baillie and Coddington [7] for the 2D Ising model on lattices up 
to L = 512. Our results are the following: 

1. The Li-Sokal bound is satisfied on the AT self-dual curve. (Indeed, for the 
direct algorithm we are able to prove this rigorously.) 

2. The bound is rather close to being sharp for a generic point on the AT self-dual 
curve. Using power-law fits for the quantity T int ,£/CH } we obtain estimates of 

2 By contrast, Wiseman and Domany [24] studied the single- cluster ("Wolff") version of the 
direct algorithm (in the single-cluster context the direct and embedding algorithms turn out to be 
equivalent). It is important that both the multi-cluster and single-cluster versions be studied, as 
they may well lie in different dynamic universality classes. We concentrate here on the multi-cluster 
version, because it is only for this version of cluster algorithms that the Li-Sokal bound (1.1)/(1.2) 
is known. 
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Zint,£ — Oi/v ranging from ss0.05 to ss0.12. The lower value corresponds to the 
Ising and X2 cases, and the higher value to the 4-state Potts model. 

3. In all cases the data are consistent with a logarithmic growth of t^^/Ch as 
A + B log L. For the Ising and the ZF models, a logarithmic behavior A log p L 
with p 0.31 also gives a reasonable fit. 

4. In all cases, the data are consistent with the boundedness of t^^/Ch as 
L — > oo only if one assumes rather strong corrections to scaling, i.e., A-\-BL~ A 
with | < A < |. Moreover, for the 4-state Potts model this scenario implies 
an implausibly large value for the coefficient B (\B\ > 10). In all cases the 
A + BL~ A fit is inferior to the AL' P and A + B log L fits. 

5. If we believe (on theoretical grounds) that there should be some continuity 
in the behavior of the ratio t^^/Ch along the self-dual curve of the AT 
model, then the possible scenarios reduce to the pure power-law and simple 
logarithmic A + B log L behaviors. 

Thus, the bound (1.1) on the autocorrelation time is not sharp, but it might be sharp 
modulo a logarithm at some or all points of the AT self-dual curve. Further studies 
on significantly larger lattices will be required in order to distinguish convincingly 
between a logarithmic and a small-power-law growth. 

This paper is organized as follows: Section 2 reviews the definition and properties 
of the AT model: phase diagram, critical exponents, etc. In Section 3 we construct 
our SW-type algorithms for the AT model, and we relate them to other SW-type 
algorithms. In Section 4 we present and analyze our numerical results for three 
selected points on the AT self-dual curve. The sharpness of the Li-Sokal bound is 
also discussed. Finally, in Section 5 we summarize our conclusions. In Appendix A 
we provide a rigorous proof of the Li-Sokal bound for the direct AT algorithm. In 
Appendix B we explain in detail how we performed the fits of the autocorrelation 
functions to extract estimates of the exponential autocorrelation times. 

2 The Ashkin-Teller model 

The Ashkin-Teller (AT) model [23] is a generalization of the Ising model to a 
four-state model. To each lattice site x we assign two Ising spins a x = ±1 and 
t x = ±1, and they interact through the Hamiltonian 

if AT = —J^(7 x (Ty — j'^2 l T x Ty — K^2 l (T x T x (T y Ty, (2.1) 

(xy) (xy) (xy) 

where the sums run over nearest-neighbor pairs (xy). It can be interpreted as two 
Ising models with nearest-neighbor couplings J and J' and interacting via a four- 
spin coupling K. Note that the fields cr, r and gt play symmetric roles in this 
model; we can consider any two of these three as the "fundamental fields". 

This model contains as particular cases some other well-known systems. The 
plane K = in the coupling-constant space (J, J', A') corresponds to a pair of 
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decoupled Ising models, one with coupling J and the other with coupling J'. At 
the other extreme, the limit K — > +00 corresponds to a single Ising model (cr = r) 
with coupling J + J'. Finally, the line J = J' = K is the 4-state Potts model with 
•/potts = 4 J . 

The family (2.1) of AT Hamiltonians exhibits several symmetries. First of all, we 
can permute freely the spin variables (cr, r, err). This implies that the AT model is 
mapped onto an essentially equivalent model under any permutation of the couplings 
(J, J',K). Moreover, if the lattice is bipartite we can flip a or r or both on one of 
the two sublattices (in other words, we choose exactly two of the three variables cr, 
r, gt and flip these on the chosen sublattice). This implies that the AT model is 
mapped onto an essentially equivalent model under the following transformations: 

(J, J', K) -> (- J, J', -K) (2.2a) 
(J, J', K) -> (J, - J', -K) (2.2b) 
(J, J', K) -> (- J, - J', A') (2.2c) 

These transformations will be useful in Section 3. 

We are mainly interested in one particular case of the Hamiltonian (2.1): the 
symmetric 3 AT model characterized by J = J', 

#SAT = — J {(7 x (T y + T x T y ) — K ^2 (T x T x (T y Ty . (2.3) 
(xy) (xy) 

This is exactly the general Z 4 clock model 4 

#dock = -2 J ]T cos (0 X - 6 y ) - K ]T cos (20 x - 26 y ) , (2.4) 

( X y) (xy) 

where the dynamical variables take the values 6 X £ {0, 7r/2, 7r, 37r/2}. The relation 
can be easily seen by setting a x = \/2cos(8 x — 7r/4) and t x = \/2cos(8 x + vr/4). 

The AT model exhibits a rich phase diagram, in both two [26,27,28] and three 
[26] dimensions. Here we concentrate on the two-dimensional square-lattice sym- 
metric AT model. Although we do not know how to solve this model analytically, 
we have a fairly good understanding of its phase diagram (see Figure 1). From 
(2.2c) we see that in the symmetric AT model a sublattice flip of a and r corre- 
sponds to the change J — > — J; it follows that the phase diagram is symmetric under 
reflection in the J = axis, under which ferromagnetic a and r ordering becomes 
antiferromagnetic (AF) and vice versa. For this reason, we show in Figure 1 only 
the half-plane J > 0. 

The line K = corresponds to a pair of decoupled Ising models, so there are 
Ising critical points at (J, K) = (±| log(l + -\/2), 0). Point DIs in Figure 1 represents 

3 Baxter [27] calls this model the "isotropic" AT model. We prefer not to use this terminology, 
in order to avoid confusion with spatial isotropy or anisotropy. 

4 More precisely, this is the general Z 4 clock model on an undirected graph. On a directed graph 
(i.e., one with oriented bonds (xy)), the interaction term sin(6 x — 6 y ) = \{(J x Ty — cr y T x ) is also 
allowed. 
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the one with the plus sign (i.e., the ferromagnetic one). The model with J = is 
again an Ising model, but in the variable err. There are thus additional critical Ising 
points (J, K) = (0, ±| log(l + \/2)). The one with the plus sign (point Is in Figure 1) 
is ferromagnetic, while the one with the minus sign (AFIs) is antiferromagnetic. 
Finally, in the limit K — > +oo we find Ising transition points at J = ±| log(l + \/2)- 
The one with the plus sign is ferromagnetic and corresponds to the point Is' of 
Figure 1. The other one is antiferromagnetic and is not depicted in Figure 1. 

The line J = K corresponds to the 4-state Potts-model subspace (right dash- 
dotted line in Figure 1). Therefore, there is a ferromagnetic critical point at J = 
K = | log 3 (point P in Figure 1). The AF regime (J = K < 0) is more subtle. 
There is no rigorous result concerning the existence or non-existence of a critical 
point in the AF 4-state Potts model. However, there is a strong numerical indication 
[29] that this model is non-critical, even at zero temperature: indeed, the second- 
moment correlation length £ is <S 2 lattice spacings at all temperatures, uniformly 
down to T = 0. 5 The absence of a critical point along the line J = —K (left 
dash-dotted line in Figure 1) follows immediately using the J — > — J invariance. 

The AT model on any planar graph can be mapped into another AT model on 
the dual graph [27,30,31]. The duality transformation is best viewed in terms of the 
Boltzmann weights 6 

io = e J+J ' +K+Jo (2.5a) 
u t = e J - J '- K+Jo (2.5b) 
u 2 = e -J+J'~K+J (2 _ 5c) 

uj 3 = e - J - J '+ K + J o ( 2 .5d) 

where J is an arbitrary constant fixing the zero of energy. The AT model with 
Boltzmann weights (cu , lo 2} 003) is mapped by duality to a new AT model with 
weights (lu , &2, £3) given by 7 ' 8 

lo = \(u + ui + u 2 + L03) (2.6a) 

ui = \(u + ui - u 2 - L03) (2.6b) 

u 2 = \{w Q - ui + u 2 - co 3 ) (2.6c) 

uj 3 = \(to Q - to x - to 2 + u 3 ) (2.6d) 



5 The critical properties of the antiferromagnetic (/-state Potts model depend strongly on the lat- 
tice structure. For instance, the AF 3-state Potts model has a transition at non-zero temperature on 
the triangular lattice [32], has a critical point at zero temperature on the square lattice [29,33,34,35], 
and is expected to be non-critical at all temperatures on the hexagonal lattice. 

6 The four energy states on a bond (xy) are labeled 0,1,2,3 as follows: the high-order bit is 
(1 — a x a y )/2 and the low-order bit is (1 — T x T y )/2. 

7 Note that if the original weights w; are normalized so that Jo = 0, the dual weights w; do not 
necessarily have this normalization. 

8 This duality transformation corresponds to the Fourier transform on Z2 x Z2 followed by inter- 
change of a and r (i.e., interchange of u>\ and u>2). 
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The symmetric AT model (the one with u\ = cu 2 ) is clearly mapped under duality 
into another symmetric AT model (i.e., uj\ = uj 2 )- Specializing to the square lattice, 
the dual graph is again a square lattice, and the self-dual manifold of (2.6) is 

coo = ui + u 2 + Lo 3 ■ (2.7) 

For the symmetric AT model on the square lattice, the self-duality condition (2.7) 
can be easily written in terms of the coupling constants: 

e~ 2K = sinh2J . (2.8) 

This is represented in Figure 1 by the curve B DIs P C. 

Furthermore, the AT model on any planar graph can be mapped onto an 8- 
vertex model on the medial graph [31]. In particular, the AT model on the square 
lattice can be mapped [27] onto a staggered 8-vertex model on the square lattice 
(which has not been exactly solved in general). As a special case, the AT model 
on the self-dual manifold (2.7) maps onto a homogeneous 8-vertex model, which is 
exactly soluble. Furthermore, the symmetric self-dual AT model (2.8) maps (after 
a simple further transformation) onto a homogeneous 6-vertex model. In this way, 
Baxter showed that the self-dual curve (2.8) is critical only for K < | log 3 (solid 
curve in Figure 1), and is non-critical for K > |log3 (dotted curve in Figure 1). 
The critical part belongs to the universality classes of the conformal field theories 
with central charge c = 1 (i.e., it can be related to the Gaussian model [36]). Along 
this line the critical exponents vary continuously and they are known exactly (see 
below). 

From series expansions [26], mean-field theory and approximate real-space renorm- 
alization-group calculations [28] , we know that two critical curves emerge from the 
Potts point P, one going to the Ising critical point Is and the other one going to the 
Ising critical point Is' at K = +00. 9 The critical curves P-Is and P-Is' map into 
one another under the duality relation (2.6) [27]. Finally, there is another critical 
curve emerging from the point AFIs and pointing towards K — > —00. The exact 
location of these three curves is unknown, as is their universality class. However, 
most people believe that they are Ising-like. In [38] it is argued that these critical 
curves should be given by non-algebraic functions. 

The four critical curves mentioned above are the borderlines of the four phases 
appearing in this model for J > 0. These phases are 10 : 

I. This is the so-called Baxter phase [26]. The spins a and r are independently 
ferromagnetically ordered. There are thus four extremal infinite- volume Gibbs 
measures according to the signs of (a) and (r) (which may be chosen inde- 
pendently); the sign of (err) is then equal to that of (cr)(r). 

9 Pfister [37] has proven the existence of two phase transitions (i.e., of the three phases II, III and I 
in succession) along any ray in the quadrant J, K > with slope < J/K < More generally, this 
applies in the full AT model along any ray in the octant J, J', K > with slope < (J + J')/K < 1. 

10 We follow the terminology used in Baxter's book [27]. 
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II. This is the paramagnetic phase, in which all three spins cr, r and gt are 
disordered. There is a unique infinite- volume Gibbs measure. 

III. In this phase both the spins a and r are disordered (i.e., lim\ x _ y ^ 00 (a x a y ) = 
YiTD.\ x _y\^ 00 {T x T y ) = 0), but their product gt is ferromagnetically ordered (i.e., 
lim^-yi^oo {cr x T x a y Ty) > 0). There are two extremal infinite- volume Gibbs 
measures according to the sign of (err). 

IV. This is the antiferromagnetic analogue of Phase III: the spins a and r are dis- 
ordered, while the product gt is antiferromagnetically ordered (i.e., lim\ x _ y ^ 00 
( — l)!^ - ^ ((7 x T x (7 y Ty) > 0). There are again two extremal infinite- volume Gibbs 
measures, according to the sign of the sublattice magnetization ( — 1)^{ct x t x ). 

This picture has been proven rigorously for low temperature and arbitrary spatial 
dimension, using Pigorov-Sinai theory [39]. In particular, it was shown that deep 
in region I there exist four periodic extremal Gibbs measures, and deep in regions 
III and IV there exist two periodic extremal Gibbs measures. 11 

The critical exponents along the self-dual curve can be obtained by relating 
the AT model to the 8-vertex model or to the Gaussian model [36,41,42]. We 
parametrize the critical part of the self-dual curve by 

e* J = ^ 2 + 2c ° S/i + 1 (2.9a) 
V2 + 2 cos \i - 1 v ' 

e 4K = l + 2cos/i (2.9b) 

where < /i < 2tt/3. This parameter /i is related to the coupling constant g of 
the Gaussian model [43] 12 by /i = 7r(l — y/4), so that | < g < 4. Thus, g = | 
corresponds to the point at K = — oo (B in Figure 1), g = 2 is the decoupled Ising 
model, g = 3 is the model considered by Zamolodchikov and Fateev [45], and g = 4 
is the 4-state Potts model. The critical exponents along the self-dual curve are 
given by 



2- y 

v = 

3- 2y 

a _ 2 - 2y 
u ~ 2-y 

1 = 1 
v 4 

Y = 7-4y 
v 4 - 2y 



(2.10a) 
(2.10b) 
(2.10c) 
(2.10d) 



11 Note the order of adjectives: Pigorov-Sinai theory studies extremal Gibbs measures that happen 
to be periodic; it says nothing about nonperiodic Gibbs measures (e.g., those with interfaces) or 
about periodic Gibbs measures that are extremal only within the restricted class of periodic Gibbs 
measures. For further discussion, see [40, Section B.3.1]. 

12 Our g corresponds to that of Saleur [43], and equals 2ir times the K of Kadanoff and Brown [36] 
and Yang [44]. 
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where the parameter y is related to /i by 

y = — = 2-f (2.11) 

with < y < |. Here a is the specific-heat exponent, while 7 (resp. 7') is the 
susceptibility exponent for a and r (resp. for err). We have chosen to specify the 
ratios (2.10b-d) because these are directly measurable by our Monte Carlo methods. 
For y > 1 (corresponding to K < 0), the specific heat has a cusp singularity (a < 0) 
rather than a divergence, but both susceptibilities remain divergent. 

The region between the decoupled Ising model (DIs) and the critical 4-state 
Potts model (P) is the most interesting. It is worth mentioning that both the 
g-state Potts model at criticality and the symmetric AT on the self-dual curve 
can be represented as certain 6-vertex models [27,41,46]. By relating these two 
6-vertex models, we can map the former model onto the latter one, and use q as 
a parametrization of this subset of the AT self-dual curve. For the square lattice, 
the g-state Potts model at criticality is mapped to the point given by (2.9) with 
2 cos /i = ^Jq. Thus, the case q = is mapped to K = (i.e., the decoupled Ising 
model), q = 2 to the model considered by Zamolodchikov and Fateev [45], and q = 4 
to the point J = K = |log3 (i.e., the 4-state Potts model). 



3 The Algorithm 
3.1 Direct algorithm 

The idea behind this new algorithm is the same as that of all Swendsen- Wang- 
type algorithms [47] 13 : we decompose the Boltzmann weight by introducing new 
dynamical variables (living on the bonds of the lattice), and we then simulate the 
joint model of old and new variables by alternately updating one set of variables 
conditional on the other set. As we have two distinct sets of Ising spins, we expect 
to introduce two distinct sets of auxiliary variables. 

We begin by enumerating the possible energy values which can occur on a given 
bond (xy). Out of the 16 spin configurations on each bond, there are only four 
different energy values (see Table 1). We can order these energies in increasing 
order, but this ordering of course depends on the relative values of the coupling 
constants J, J' and K . Instead of developing a different algorithm for each possible 
ordering, we can use the symmetries (2.2a-c) of the general AT Hamiltonian (2.1) 
and choose 14 an equivalent general AT model satisfying 

J,J'>\K\. (3.1) 

For such a model, the energies satisfy 

= E < E U E 2 < E 3 . (3.2) 

13 For a pedagogical presentation, see [1, Section 6] or [4, Section 4]. 
14 At least if the lattice is bipartite. 
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In particular, for a given bond, the lowest-energy state is the one with the two a 
spins parallel and the two r spins parallel. We shall hereafter assume that (3.1) 
holds. 

Remark. The following algorithm is also valid for a non-homogeneous AT model 
on an arbitrary finite graph, with Hamiltonian 

Hat = —^2 JxyVxVy — ^2 J'xy T x T y ~ ^2 K xy a x T x a y T y (3.3) 

(xy) (xy) (xy) 

satisfying 

Jxy,J xy ^ \^xy\ (3-4) 

for every bond (xy). It suffices to make the obvious notational alterations. 

The Boltzmann weight associated with the bond (xy) is equal to 

W hoil d(cr x ,ay,T x ,Ty) = e " 2(J+J,) + 

e~ 2J ' [e~ 2K ~ e" 2J ] K x , y + 
e~ 2J [e~ 2K ~ e- 2J 'j S Tx , Ty + 

[l " e-« J ' +K ) - e-W) + e- 2 ( J+J ')] <W Wa . (3.5) 

Let us now introduce the auxiliary variables m xy} n xy = 0,1 associated with the 
spins cr, t respectively, and define the joint-model Boltzmann weight for the bond 
(xy) to be 

Wi°^((T x ,(Ty,T x ,Ty;m X y,n X y) = e - 2 ( J+ J ' ) 6 m x y , o 6 nxy , o + 

-23' \ -2K -23] c c c , 

e [e - e j c ax , ay b mxyA b nxyfl + 
e" 2J [e~ 2K ~ e~ 2J '\ S Tx)Ty 6 mxy)0 6 nxy)1 + 

[l " e-W) - e~ 2 ^ + e- 2 ^} 6^ y 6 Tx , y 6 mxyA 6 nxyA . (3.6) 

If we sum (3.6) over m xy and n xy , we obtain (3.5), as desired. The complete joint 
probability is then 

W^({a } r}; {m, n}) = j Wg£(v x , a y , r X} r y ; m xyi n xy ) , (3.7) 

xy 

where Z is the partition function of both the joint model and the original model: 
Z = Ii W bZd( (J x,(Ty,T x ,Ty;m X y,n X y) (3.8a) 

<T,T = ±1 771,71 = 0,1 (xy) 

= Y[W hond (a x ,a y ,T x ,Ty) (3.8b) 

<T,T = ±1 ( X y) 

Our SW-type algorithm consists in simulating the joint probability distribution 
(3.7) by alternately applying the conditional distributions on {cr, r} and on {m,n}. 
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These two steps can be read off immediately from (3.6)/(3.7); in detail, they are 
the following: 

Step 1: Update of {m,n} given {cr, r}. Conditional on the {cr, r} config- 
uration, the bond variables {m, n} are given independently for each bond. For a 
bond (xy) with spins g X} a y} t X} t V7 we obtain the new bond variables m xy and n xy 
(independently of the old values) by the following rules: 

la) If o x = G y and t x = T y} then we choose (m xy} n xy ) with the following proba- 
bilities 

(m xy ,n xy ) = (1, 1) with pi = 1 - e- 2 ( J '+ A ") - e " 2 ( J +^) + e - 2 (^') 



(m xy ,n xy ) = (1,0) with p 2 
(m^, n xy ) = (0, 1) with p 3 



2J' L-2A' „-2J 



e e — e 



2J „-2A' „-2J' 



e e — e 



(m xy , n xy ) = (0, 0) with p 4 = e 2 ( J+J ') = 1 - Pl - p 2 - p 3 . 

lb) If a x = G y and t x = —T y} then the probabilities are 

(m X y,n X y) = (1,0) with q 1 = 1 - e~ 2 ( J+K K 
(m X y,n xy ) = (0,0) with q 2 = e -^ J+K "> = 1 - q 1 . 

lc) If g x = — G y and t x = T y} the probabilities are 

(m xy ,n xy ) = (0, 1) with n = 1 - e- 2 ( J '+ A "). 
(m xy ,n xy ) = (0,0) with r 2 = e - 2 ( J ' +A ") = 1 - r x . 

Id) If g x = —G y and t x = —T y} we choose (m xy} n xy ) = (0,0) with probability 1. 

All these choices are made independently for each bond (xy). 

Step 2: Update of {cr, r} given {m, n}. Given the bond-configuration {m, n}, 
we build all the connected clusters of g spins (resp. r spins) joined by bonds with 
m xy = 1 (resp. n xy = 1). Within each cluster, the spin values are required to be 
equal, but this common value may be either +1 or —1. The spin value for each 
cluster is chosen randomly, independently of the old value and of the choices made 
for the other clusters. 

One iteration of the direct algorithm consists of an application of Step 1 followed 
by an application of Step 2. 

Remarks. 1. Wiseman and Domany [24] have introduced essentially this same 
decomposition of the Boltzmann weight, although their derivation is in our opin- 
ion more complicated. They then studied numerically the single-cluster ("Wolff") 
version of this algorithm. Here we study the many-cluster ("Swendsen-Wang") 
version. 

2. This direct SW-type algorithm satisfies the Li-Sokal bound (1.1)/ (1.2). The 
proof is a straightforward generalization of the one given in [9] for the Potts case; 
we present it in Appendix A. 
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3. We can generalize our SW-type algorithm to a "generalized Ashkin-Teller 
model" consisting of a g-state Potts variable a and an r-state Potts variable r 
interacting through the Hamiltonian 

H GAT = -2( J - K) J2 " 2( J' - K) £ 8 t ,t v ~ ±K £ S^ y S TxTy . (3.9) 

(xy) (xy) (xy) 

It is clear that from this Hamiltonian we obtain again the joint probability distri- 
bution (3.6). This model has been considered in [25] when J' = K . 

3.2 The embedding algorithm 

The algorithm presented in the preceding subsection is perfectly legal, but it is 
somewhat complicated to write the computer code for its Step 1 in an efficient way. 
In this section we introduce a variant algorithm in which we deal with only one 
kind of spin (a or r) at a time. 

Consider the Boltzmann weight of a given bond (xy) } conditional on the {r} 
configuration (i.e., the r spins are kept fixed): it is 

W hond (a X} a yiTxi Ty ) = e" + [l - e " W^,)] S<7xj<7y . (3.10) 

We can simulate this system of a spins using a standard SW algorithm. The effective 
nearest-neighbor coupling 

Jf y = J + Kr x r y (3.11) 

is no longer translation-invariant, but this does not matter. The key point is that 
the effective coupling is always ferromagnetic, due to the condition (3.1). An exactly 
analogous argument applies to the {r} spins when the {a} spins are held fixed. 
The embedding algorithm for the AT model has therefore two parts: 

Step 1: Update of {a} spins. Given the {r} configuration (which we hold 
fixed), we perform a standard SW iteration on the a spins. The probability p xy 
arising in the SW algorithm takes the value p xy = 1 — exp[— 2(J + Kt x t v )]. 

Step 2: Update of {r} spins. Given the {a} configuration (which we hold 
fixed), we perform a standard SW iteration on the r spins. The probability p xy 
arising in the SW algorithm takes the value p xy = 1 — exp[— 2(J' + Ka x a y )]. 

One iteration of the embedding algorithm consists, by definition, of a single 
application of Step 1 followed by a single application of Step 2. 

Wiseman and Domany [24] also constructed a embedding version of their single- 
cluster algorithm. Furthermore, they showed that, in the single-cluster context, the 
direct and embedding algorithms define the same dynamics 15 ; only the computer 
implementation is different. However, this equivalence does not hold for our many- 
cluster algorithm. In the direct algorithm we have independent clusters of a spins 

15 More precisely, this equivalence holds when the embedding algorithm is defined by making a 
random choice of Step 1 or Step 2 at each iteration. 
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and t spins that could be flipped simultaneously. In the embedding algorithm we 
have at each step only one of the two types of clusters. 

The embedding algorithm, due to its simplicity, is the one used in our MC study 
of the AT model (see Section 4). 

Remark: This algorithm is closely related to an embedding algorithm for general 
Z n clock models on an undirected graph. We can consider Z n as a subgroup of U(l) 
and then apply Wolff's embedding algorithm for the XY model [12,48,49]. Let us 
specify the reflection plane by a vector s = (cos </>, sin cf>) in this plane; clearly <f> 
is specified only modulo tt. If n is odd, there is a unique type of reflection: the 
reflection plane passes through one spin value and one point bisecting two spin 
values, and it corresponds to <f> = 2Trk/n with k either integer or half-integer. How- 
ever, if n is even, there are two types of reflections: the reflection plane can either 
pass through two spin values or else through two bisector values; these correspond 
to <f> = 2Trk/n with k integer or half-integer, respectively. Thus, for the 4-state 
clock model we have two reflections of the first type and two of the second type: 
with the identifications a x = \/2cos(8 x — 7r/4) and t x = \/2cos(8 x + 7r/4) [taking 
6 X £ {0, 7r/2, 7r, 37r/2}], these reflections are 

^ = 0: (a, r) -> (r, a) . 

<P = tt/2: (a,T)^(-T,-a) ^ AZ > 

and 

^ = tt/4: (<7, r ) ^(<r,-r) V ' ; 

respectively. The last two moves (i.e., those of bisector type) are precisely the moves 
allowed in our algorithm. The first two moves correspond to the interchange of a 
and r, either without or with a simultaneous flip of both spins. 

So let us fix <f> to be one of the four values listed above, and let us embed Ising 
spins e x = ±1 into the 4-state clock model via the Wolff update 

e x $ + s x (6 x - <)>) . (3.14) 

Plugging this into the clock-model Hamiltonian (2.4), we obtain an induced Ising 
system in the {e} variables with interaction 

Jf y = 2 | sin(^ - <j>)sm{e y - <f>)\ [J + 2K cos{6 x - <t>)cos(6 y - <f>)] (3.15) 

This Ising system can then be simulated by means of an SW algorithm. The cases 
<f> = T 71 "/ 4 correspond to Steps 1 and 2 of our embedding algorithm. On the other 
hand, the case <f> = is characterized by an effective coupling 

Jf y = 2JS^_ T J^_ Ty . (3.16) 

This means that only bonds joining sites with a x ^ t x and a y ^ r y can be "ac- 
tivated" [with probability p = 1 — exp(— 4J)]; the move (cr, r) — > (r, a) is then 
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equivalent to flipping both a and r within each such cluster. An analogous conclu- 
sion applies to the case <f> = tt/2: here only bonds joining sites with a x = t x and 
G y = T y can be "activated". Thus, the moves with <f> = 0,7r/2 are in a sense merely 
combinations of the moves <f> = ±7r/4 already contained in our AT embedding al- 
gorithm. For this reason, we think that the introduction of the moves <f> = 0,7r/2 
into our embedding algorithm will not further reduce the dynamic critical exponent. 
Indeed, the algorithm with only the <f> = 0,7r/2 moves is not even ergodic: at each 
site the product (t x t x is conserved. 

3.3 Particular cases 

As pointed out in Section 2, the AT model reduces to two decoupled Ising models 
at K = and to the 4-state Potts model al J = J' = K . It is worth mentioning 
that the above-discussed algorithms reduce to the well-known SW algorithms for 
those particular cases. 

When K = 0, it is easy to verify that (3.6) reduces to 

WbondC* 7 ** a v, T ^ r v\ ™>xvi n xy) = [e~ 2J 6m xy ,o + (l - e" 2J ) <W a <Wi] x 

[e" 2J '^,o + (l - e- 2J ') <Ww] • (3-17) 

This means that the Boltzmann weight for any bond is just the product of the 
weights of the two independent Ising models. As a result, our direct AT algorithm 
reduces to two independent SW algorithms on the systems {cr, m} and {r, n}. Of 
course, the same holds for the embedding algorithm, as the a spins are decoupled 
from the r spins. 

When J = J' = K > 0, (3.6) can be written as 

^bo'ncK ^, &yi T x , T y ] m xy , n xy ) = e 4J fim xy ,ofin xy ,0 + 

[! - e " 4j ] S<7 X ,<7yST X ,TyS mX y,lSn X y,l , (3-18) 

which is exactly the standard SW decomposition for the Boltzmann weight of the 
4-state ferromagnetic Potts model. As a result, our direct AT algorithm on the line 
J = J' = K > reduces to the standard SW algorithm for the 4-state ferromagnetic 
Potts model. However, the embedding algorithm does not reduce to the standard 
SW algorithm in this case. 

4 Numerical Results 

4.1 Autocorrelation functions and autocorrelation times 

We are interested in the dynamic behavior of the embedding SW algorithm 
described in Section 3.2. Thus, we need to study the autocorrelation functions and 
autocorrelation times for each measured observable. Given an observable 0, we 
define the corresponding unnormalized autocorrelation function as 

Coo(t) = (O s O s+t ) - {Of , (4.1) 
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where all the expectation values (•} are taken in equilibrium and t is the "time" in 
units of MC steps. 16 The associated normalized autocorrelation function is 

Poo(t) = Coo(t)/Coo(0) . (4.2) 
The integrated autocorrelation time for the observable O is defined as 

2 CO 

Tint,© = - POo(t) 
Z t=-oo 
2 oo 

= li + T.Pooit). (4.3) 

z t=i 

Here the factor | is purely a matter of convention (if the normalized autocorrela- 
tion function is a pure exponential, poo{t) ~ e~'*'/ T with r ^> 1, then this definition 
implies that T int ,o ~ t). Finally, the exponential autocorrelation time for the ob- 
servable O is defined as 

T exp ,e> = lim sup — ^—jy- — , (4.4) 



t^OO 



log \poo(t)\ 



and the exponential autocorrelation time ("slowest mode") for the system as a whole 
is defined as 

r exp = sup Te XPi e> • (4.5) 
o 

Note that r exp = T exPj £> whenever the observable O is not orthogonal to the slowest 
mode of the system. 

The integrated autocorrelation time controls the statistical error in Monte Carlo 
estimates of the mean (O). In particular, given a sequence of n Monte Carlo mea- 
surements of the observable O — call them {0\ } . . . , O n } — the sample mean 

— 1 n 

O = - Y,Ot (4.6) 
n 7^i 



h 



as a variance 



var 



1 n 

(O) = - E C 'oo(r ~ s) (4.7a) 



n 2 1 

r,s=l 



= ~ E ^--)coo(t) (4.7b) 

~ - 2r intj o Coo(0) for n > r intj o (4.7c) 
n 

This means that the variance is a factor 2T int ,o larger than it would be if the mea- 
surements were uncorrelated. It is therefore, very important to estimate the au- 
tocorrelation times for all the interesting observables in order to ensure a correct 



16 0ne "Monte Carlo step" consists of one application of "step 1" of Section 3.2 followed by one 
application of "step 2" . 
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determination of the statistical errors. The integrated autocorrelation time T int ,o 
can be estimated using standard procedures of statistical time-series analysis [50,51] . 
In this way, we obtain reliable error bars for both T int ,o and (O). We have used a 
self-consistent truncation window of width 6T int ,o [52, Appendix C]. This window 
width is sufficient whenever the autocorrelation function pooif) decays roughly ex- 
ponentially, a behavior that we will confirm explicitly here (see Section 4.5). 

The exponential autocorrelation times T exPj e> are extracted by fitting the autocor- 
relation function pooif), f° r t large enough, to a pure exponential A exp(— t/T exPj e>). 
The statistical details of this fit are described in Appendix B. However, it should 
be emphasized that it is in principle impossible to obtain a statistically valid esti- 
mate of T exPj £>, as there could always be a very- slowly- decaying component of pooif) 
with arbitrarily small amplitude (which would thus be invisible under the statistical 
noise). Thus, our estimates of T exPj e> are really lower bounds. They can be taken as 
estimates of T exPj e> only if one assumes that pooif) 1S roughly an exponential, with 
no very- slowly- decaying components. 



4.2 Observables to be measured 



Let us begin by defining some basic observables. The observables of interest 
involving only the a spins are 



E 



X 



a me 



2ttx 1 /L 



+ 



X 



cr me 



2ttx 2 /L 



(4.8) 
(4.9) 

(4.10) 



where L is the linear size of the system (we always use periodic boundary condi- 
tions) and (xi, x 2 ) are the Cartesian coordinates of the point x. The last observable 
can be also seen as the square of the Fourier transform of a at the smallest al- 
lowed non-zero momenta (i.e., (±27r/L,0) and (0,±27r/L) for the square lattice); 
it is normalized to be comparable to its zero-momentum analogue M. a . We define 
analogous observables for the r spins and for the composite operator err. 

We have run the MC algorithm on three different points of the self-dual curve of 
the symmetric AT model (see Table 2). One is the 4-state Potts model at criticality, 
where the three variables (cr, r, gt) are related by symmetry. But for the rest of 
the points of the self-dual curve, only a and r are related by symmetry. Since we 
wish to exploit the symmetries of the model in our data analysis, our choice of 
observables to measure will depend on which model we are studying. 



4.2.1 Observables for the critical 4-state Potts model 

For AT models on the 4-state Potts line J = J' = K , the natural choice of the 
observables are those having the symmetries of the original Potts model, namely 
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those invariant under permutations of (cr, r, err). We have measured the expectations 
and autocorrelation times for the following observables: 



M 2 = l(Ml + M 2 T + Ml T ) (4.11) 

£ EE | (E a + £ T + Ear) (4.12) 

T = \(T a +T T + T aT ) (4.13) 

These observables coincide with the usual ones for the 4-state Potts model up to 
some multiplicative constants. We can then define the magnetic susceptibility 

X = \{M 2 ) , (4.14) 



the total energy 17 
the specific heat 
and the second-moment correlation length 



E = ^(S) , (4.15) 



Ch = {(£ 2 ) - (£) 2 ) , (4.16) 



where F is defined as 



F . I (4.17) 

2 sin f K J 



F = i(^) . (4.18) 



In all these formulae, V is the number of lattice sites (i.e., V = L 2 ) and 2V is 
the number of bonds. This definition of the correlation length is not equal to 
the exponential correlation length (= 1/mass gap), but it is expected that both 
correlation lengths scale in the same way as we approach the critical point. 

Remark. To compute the error bar of the specific heat Ch we have first computed 
the mean energy (£}, and then considered the observable O = {£ — (S)) 2 using the 
procedures described in this section. 



4.2.2 Observables for a general point on the self-dual curve 

A generic point on the self-dual curve does not enjoy the 4-state Potts symmetry, 
but it does of course enjoy the a <->■ r symmetry characteristic of all symmetric AT 
models. The natural choice is thus to define two sets of observables: one for the a 

17 We have normalized the energy such that —1<E<1 (i.e., E is the energy density per bond), 
and the same normalization has been taken for the specific heat. However, in the literature it is 
more common to find the energy and specific heat normalized per site, i.e., with a factor l/V rather 
than our l/(2V). 



18 



and t variables, and another for the composite operator err. The first set is given 
by 



Ml = \{Ml + Ml) (4.19) 
£ w = \(£ a + £ T ) (4.20) 
T w = \{? a + ? T ) (4.21) 



and the second one by M aT , £ aT and T aT . We then define 



- l 1/2 

X.= y(M 2 J, E U = -^(E U ), F U = ±(F U ), (4.22) 



and 



1 , . .o . „ 1 ,„ . „ 1 



1/2 



X(jt — ^.(A / ( CTT ) , Ear — 2y(^' <7T ) ' -^ ,<7T — y(*^ 7<7T ) ' l ? <7T ~~ 2 sin — 
Finally, we define a specific-heat matrix Ch 



(4.23) 



Pi _ J_ ( var(£ w ) cov(^,£ (7T ) \ . . 

H ~ 2V\cov(£^£ aT ) var^) J ' 

with eigenvalues C£f jin i n and Ch, max an d corresponding eigenvectors 

( i "\ 



^min I 

\ a 



(4.25a) 

^max = 1^1 (4.25b) 



a 



where a is some real number. These two eigenvalues have distinct critical exponents: 
Cff,max corresponds to a relevant operator (with exponent a > 0), while CH,min is 
expected to correspond to a marginal operator (with exponent 0). The marginal 
operator arises from the existence of the self-dual curve (2.8), along which the 
critical exponents vary continuously. In particular, on the self-dual curve (2.8), we 
expect that 

a = -±coth2J (4.26) 

(in the infinite- volume limit), as can easily be computed from the tangent vector to 
(2.8), taking into account the normalization £ = 2J£ LU + K£ aT . 



4.3 Summary of our simulations 

We have run our MC program for the embedding algorithm (Section 3.2) at 
three different points of the self-dual curve (2.8): see Table 2 for details. One of the 
points is the critical point of the 4-state Potts model (i.e., J = J' = K = |log3). 
The second point is the image of the q = 2 Potts model via the transformation 
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discussed at the end of Section 2 [i.e., (2.9a-b) with 2 cos /i = y/2]: this point will 
be denoted as ZF and corresponds to the model studied by Zamolodchikov and 
Fateev [45]. The third point is one of the ones studied in [24] and will be denoted 
as in that paper (X2); it corresponds to a Potts model with q 0.651287. We also 
notice that the point X3 of [24] is rather close to our choice ZF. Finally, we have 
used the extensive MC data of Baillie and Coddington [7] 18 for the critical Ising 
model (which corresponds to the point DIs of the AT model). In Table 3 we include 
the static and dynamic data corresponding to this point. 

For each of these points we have run the MC program at different lattice sizes 
ranging from L = 16 to L = 1024 for the 4-state Potts model and from L = 16 
to L = 512 for the other two points. In all cases we have started the simulations 
with a random configuration and we have discarded the first 10 5 iterations to allow 
the system to reach thermodynamic equilibrium. This discard interval is sufficient 
for equilibration: in the worst case (i.e., the 4-state Potts model with L = 1024), 
it is roughly equal to 190ri ntj £ (or about 160r exPj £). The number of measurements 
ranges between 8 X 10 5 and 4.4 X 10 6 . In all cases except the L = 1024 Potts, the 
number of measurements is greater than 10 4 Ti ntj £ . This is sufficient to obtain good 
estimates (errors ~ 1-4%) for the autocorrelation times, and excellent estimates 
(errors ~ 0.1-1%) for the static observables. On the other hand, for the 4-state 
Potts model with L = 1024 we were able to achieve only ~ 1500ri ntj £. The error 
bars on this point are therefore rather large. 

To test the program we have compared the MC results to the exact solution 
for small lattices (3x3 and 4x4). We performed this test over a wide range of 
couplings (J, J', A'), including both high- and low-temperature regions as well as 
the critical region. We have also compared the results for larger lattices to previous 
MC computations [9] for the critical 4-state Potts model. In all cases the agreement 
was good. 

The CPU time required by our program is approximately 10L 2 /is/iteration on 
an IBM RS-6000/370. The total CPU time used in this project was approximately 
1.2 years on this same machine. 

The estimates for the observables discussed in Section 4.2 are shown in Tables 4- 
8. In all of them, the quoted errors correspond to one standard deviation (i.e., 
confidence level of 68%). In reporting the number of measurements (MCS), we 
have already subtracted the MC steps discarded for equilibration. 

A summary of our estimates of critical exponents, together with the theoretically 
predicted exponents, can be found in Table 9. The data analysis leading to these 
estimates is the subject of the next two subsections. 

18 We thank Paul Coddington for communicating to us the numerical values of these data, which 
formed the basis for the graphs in [7]. For the lattices L < 128 these data coincide with the data 
reported by Baillie and Coddington in an earlier paper [6], while for L = 256, 512 they improve the 
statistics somewhat. 
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4.4 Static quantities 



For each quantity (9, we carry out a fit to the power-law Ansatz O = AL' P using 
the standard weighted least-squares method. As a precaution against corrections to 
scaling, we impose a lower cutoff L > L m i n on the data points admitted in the fit, 
and we study systematically the effects of varying L m i n . In general, our preferred 
fit corresponds to the smallest L m i n for which the goodness of fit is reasonable (e.g., 
the confidence level 19 is > 10-20%), and for which subsequent increases in L m i n do 
not cause the % 2 to drop vastly more than one unit per degree of freedom. 



4.4.1 Susceptibilities 

We have fitted the values of the susceptibilities to the power-law functions Xu> = 
AL 1 ^ and x <7T = AL 1 l v as described above. The estimates for 7/1/ and 7'/^ are, 
in all cases, very stable as we vary L m i n . This means that the corrections to scaling 
for these observables are not statistically significant (to the degree of accuracy we 
have attained here). 

For the 4-state Potts point (Table 4), our preferred estimate is obtained for 

T . _ 1«.20 

1(P) = ^(P) = 1.744 ± 0.001 (4.27) 

with x 2 = 2.21 for 5 degrees of freedom (DF), confidence level = 82%. The difference 
from the exact result (j/u = |) is small, but it is roughly equal to 6 standard 
deviations. Possibly this is due to a small correction-to-scaling effect (which has 
become statistically significant due to the very high precision we have obtained on 
the smaller lattices). For L m i n = 128 we have 

1(P) = ^(P) = 1.744 ± 0.004 (4.28) 

with x 2 = 1-54 (2 DF, level = 46%), which is now fully compatible with the exact 
result. 

For the ZF point (Table 5), our preferred estimates are obtained for L m i n = 128: 

-(ZF) = 1.750 ± 0.004 (4.29) 
v 

with x 2 = 1-54 (1 DF, level = 22%), and 

— (ZF) = 1.668 ±0.005 (4.30) 

v 



19 "Confidence level" is the probability that \ 2 would exceed the observed value, assuming that 
the underlying statistical model is correct. An unusually low confidence level (e.g., less than 5%) 
thus suggests that the underlying statistical model is incorrect — the most likely cause of which 
would be corrections to scaling. 

20 We will hereafter use the indices P, ZF and X2 to designate the points where the results apply. 
The index DIs will be used to denote the Ising model. 
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with x 2 = 1-57 (1 DF, level = 21%). The agreement with the exact results ("f/v = | 
and j'/u = §) is extremely good, and the \ 2 is acceptable. 

For the point X2 (Table 7), we get our preferred estimates for L m i n = 32: 

^(X2) = 1.750 ±0.001 (4.31) 

with x 2 = 0.98 (3 DF, level = 81%), and 

^-(X2) = 1.605 ± 0.001 (4.32) 

with x 2 = 1-24 (3 DF, level = 74%). The agreement with the exact values j/u = | 
and ^' jv = 1.6045 is again extremely good, as is the % 2 . 

Finally, we have re-analyzed the MC data of Baillie and Coddington [7] for the 
Ising model (Table 3). Our preferred fit is for L m i n = 64: 

^(DIs) = 1.7501 ± 0.0002 (4.33) 

with x 2 = 0.75 (3 DF, level = 86%). Both the high accuracy of the result and the 
agreement with the exact answer (j/u = |) are remarkable. 



4.4.2 Specific Heat 

Here we have to distinguish between the 4-state Potts model and the other 
two points (ZF and X2). For the latter points the analysis is a little bit more 
complicated, as we have to deal with a specific-heat matrix Ch [cf. (4.24)] instead 
of a single number. 

For the 4-state Potts model (Table 4) we first tried to fit the data to a pure 
power-law function Ch = AL a l v . The results of this fit (as a function of L m i n ) are 
contained in Table 10. We observe a systematic trend toward higher values of ajv 
as we increase L m i n . For L m i n > 64 the % 2 values are acceptable. Nevertheless, 
being conservative, we take as our preferred fit L m i n = 128 (boldfaced in Table 10): 

OL 

-(P) = 0.768 ± 0.009 (4.34) 
v 

with x 2 = 1-40 (2 DF, level = 50%). Clearly, the agreement with the exact known 
result {ajv = 1) leaves something to be desired! A similar result was reported by 
Wiseman and Domany [24]: ajv = 0.747 ± 0.003, using lattices 16 < L < 128. As 
a matter of fact, if we fit our own data restricted to the interval 16 < L < 128, we 
obtain ajv = 0.741 ± 0.004 ( X 2 = 3.07, 2 DF, level = 22%), which is consistent 
with the value of Wiseman and Domany. Thus, as we go to larger lattices we 
obtain estimates of ajv that are closer to the exact value, but the improvement 
from L max = 128 to L max = 1024 is extremely slow. 

However, we already know on theoretical grounds [16,17,18] that the true leading 
behavior of the specific heat involves a multiplicative logarithmic correction 

C H ~ Llog- 3/2 L. (4.35) 
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In the range of L considered here, the logarithmic factor could be mimicked by 
a power-law function, thus yielding an "effective" critical exponent (a/v)^ lower 
than ajv = 1 (in agreement with our numerical results). To check this, we tried 
to fit our data to the function Ch = AL a l v log -3 / 2 L (see Table 11). We observe 
that the estimates for ajv are much closer to the exact value. These estimates are 
slightly higher than 1, but there is a clear systematic trend towards smaller values 
of ajv as L m i n is increased. For L m i n = 128 we obtain a reasonable fit with 

OL 

-(P) = 1.039 ± 0.009 (4.36) 

with x 2 = 0.45 (2 DF, level = 80%). However, this value still differs from the true 
one by four standard deviations. 

Alternatively, we can fit the data to the trial function Ch = AL log~ p L (see 
again Table 11). We also observe a systematic trend towards the exact value p = 
| = 1.5, but the estimates of p are not compatible within errors with the exact 
value. Our preferred fit occurs for L m i n = 128, 

p(P) = 1.286 ±0.052 (4.37) 

with x 2 = 0.34 (2 DF, level = 84%). This estimate is again four standard deviations 
away from the expected result. 

In summary, it is extremely difficult is to obtain a reliable estimate of ajv when 
the leading term of the specific heat behaves like (4.35). For lattices up to L = 1024 
it is impossible to see anything even resembling the correct exponent ajv = 1. 
On the other hand, if we introduce in the fits some theoretical information (e.g., 
the power of the logarithmic term) the estimate improves a lot, but it is still four 
standard deviations away from the expected result. We need to go beyond L = 1024 
to disentangle the true asymptotic behavior of the specific heat. 

For the ZF and X2 models, we have to deal with a specific-heat matrix Ch 
as in (4.24). First we computed all its matrix elements; then we diagonalized it 
to obtain the eigenvalues CH,min and CH,ma,x and the corresponding eigenvectors 
^Vin = (1 ? °0 & nd w meix = (a, — 1) [cf. (4.25)]. The error bars on the eigenvalues 
and on the eigenvector parameter a are computed by using the standard error- 
propagation formulae. The minimum eigenvalue CH,min is expected to tend to a 
finite constant as L — > oo (i.e., to have critical exponent zero), as there should be a 
marginal operator responsible for movements along the critical self-dual curve (2.8). 
The maximum eigenvalue CH,ma,x is expected to grow as L a l v with the power given 
by (2.10b). In general the value of the parameter a varies with the lattice size L, 
and in the limit L — > oo we expect a to tend to the value 

aoo = -|coth2J (4.38) 

[cf. (4.26)]. In Table 12 we show the evolution with the lattice size of the parameter 
a corresponding to the models ZF and X2, while in Tables 5 and 7 we report the 
eigenvalues C// iinax and C// jmin . 
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ZF model. If we try to fit CH,ma,x for the point ZF (Table 5) to a pure power-law 
function, we obtain estimates of ajv which are quite consistent among themselves 
for L m i n > 64. The fit for L m i n = 64 gives 

OL 

-(ZF) = 0.663 ±0.006 (4.39) 

with x 2 = 1-12 (2 DF, level = 57%). This result agrees well with the the theoretical 
prediction ajv = |. 

The minimum eigenvalue of Ch is almost constant within statistical errors for 
L > 128. As a matter of fact, the fit to a constant C^ min is very good for L m i n = 128: 
the result is C^ min = 0.562 ± 0.004 with X 2 = 0.55 (2 DF, level = 76%). 

Finally, from Table 12 we see that the eigenvector parameter a appears to be 
tending as L — > oo to the predicted exact value (4.38). We can test this convergence 
quantitatively, and also extract an estimate of the correction-to-scaling exponent 
A, by fitting the data to a — = AL~ A . We obtain a reasonable fit already for 

A(ZF) = 0.715 ± 0.008 (4.40) 
with x 2 = 1-08 (4 DF, level = 90%). 

X2 model. The estimate of ajv for the point X2 (Table 7) is not well stabilized: 
it decreases systematically as L m i n increases (see Table 13), and for L m i n < 64 the 
X 2 values are horrible. Our preferred fit is obtained for L m i n = 128: 

OL 

-(X2) = 0.438 ± 0.008 (4.41) 
v 

with x 2 = 0.32 (1 DF, level = 57%). Notice that the % 2 value is now very reasonable. 
This estimate is still three standard deviations away from the exact result ajv = 
0.4183, but the trend is in the right direction. Moreover, for L m i n = 256 we obtain 
a slightly lower estimate, ajv = 0.430 ± 0.017, which is now consistent with the 
the exact result. The smallness of the difference between the observed and the 
true values (less than 0.02) suggests that multiplicative logarithmic corrections (as 
occur in the Potts case) are absent at this point; we are most likely seeing the 
effects of additive corrections to scaling of the form L~ A with A on the order of 
0.5 (or conceivably 1/logL). Wiseman and Domany [24] reported a value ajv = 
0.542 ±0.008, which is much larger than the exact one and than ours. This is surely 
due to the fact that they considered only rather smaller lattices (16 < L < 128). 
To check this, we fit the subset of our data corresponding to 16 < L < 128, and 
obtained a/u = 0.502 ± 0.003 with X 2 = 8.81 (2 DF, level = 1.2%). 

The smallest eigenvalue grows is again consistent with a constant C^ min within 
statistical errors for L > 128. For L mm = 128 the result is C^ min = 0.303 ± 0.002 
with x 2 = 1-17 (2 DF, level = 56%). 

Finally, the behavior of the eigenvectors is rather similar as the ZF case: they 
approach the predicted exact value (4.38) as L grows. If we fit the data to a — = 
AL~ A , we obtain for L m i n = 128 the value 

A(X2) = 0.518 ± 0.024 (4.42) 
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with x 2 = 0.02 (1 DF, level = 88%). This value for the exponent A agrees well with 
the rough estimate obtained from the corrections to scaling in the specific heat. 

Ising model. We have also computed the effective exponents (a/v)^ associated 
with a pure-power-law fit to the specific heat of the Ising model, for various intervals 
of L. Such an effective exponent will be useful as a standard of comparison for the 
numerically extracted estimates of the dynamic critical exponent z- int ,£ (see Section 
4.5.1 below). First we computed the exact values of the specific heat for a finite 
lattice, using the formulae of Ferdinand and Fisher [53]: see Table 3. We then 
performed a power-law fit, with the (fake) error bars chosen so as to give all points 
the same statistical weight; we took L max = 512 and considered various values of 
L m i n . As expected, (a/v)^ is not stable: it decreases as L m i n grows, ranging from 
(a/^) e ff = 0.244 for L m i n = 8 to (a/v)^ = 0.162 when L m i n = 256. 



4.4.3 Second-moment correlation length 

The quantity x(L) = £(L)/L is expected to approach a constant x* as L — > 
oo. This constant is characteristic of the massless Ashkin-Teller field theory on a 
continuum torus with aspect ratio 1, and in principle it should be calculable via 
conformal field theory (although to our knowledge this calculation has not yet been 
done). 

Our Monte Carlo data are consistent with this behavior. First we tried to fit 
the values for L > L m i n to a constant, using the weighted least-squares method. In 
Table 14 (fits marked C) we report our best estimates for x*. 

It is intriguing to note that the values of are all quite close to 1, though the 
differences from 1 are 7-10 standard deviations for the points ZF and X2. This 
(near-)agreement might be due to the fact that all these three models have the 
same central charge (c = 1). However, a detailed study is needed to understand the 
observed variations. On the other hand, the value of x* T manifestly decreases as we 
move towards the point K = 0, where the a spins are decoupled from the r spins. 
At that point we expect x* T = 0. 

We can also study the rate at which £(L)/L approaches x* for these three models. 
For the 4-state Potts model we already know [54] the answer for the case of the 
exponential correlation length (^ exp = 1/mass gap) on an L X oo cylinder. The 
behavior for large L is 

U P (L) 4 24 1 / 1 \ 

= ^-^t^2l + {t^l) ■ (4 - 43) 

Of course, there is no guarantee that the asymptotic behavior of Second-moment on a 
torus is the same as that of ^ exp on a cylinder, but it is a plausible guess. Therefore, 
we tried to fit our data to the function 

« £) / £ = ** + i^I • (4 - 44) 

The result can be found in Table 14 (fit marked L), and in more detail in Table 15. 
The x 2 °f this fit is excellent, but roughly the same values of % 2 could have been 
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obtained using corrections of the type 1/a/L or 1/L (see Table 15). The values 
of x* vary a little bit from one fit to the next (in some cases by several standard 
deviations). For a full understanding of these results, it would be very helpful to 
have a theoretical prediction analogous to (4.43) for the second-moment correlation 
length £ on a torus. 

For the rest of the self-dual curve we have no theoretical hint analogous to 
(4.43), so we must fit our data empirically. For the ZF model the fits of x w and 
x<j T to a constant are so perfect (L m i n = 16) that no further conclusions could be 
obtained (see Table 14). For the X2 model, by contrast, the fits to a constant need a 
much higher value of L m i n} indicating that the corrections to scaling are statistically 
significant. First we tried a correction of the type ; x = x* + AL~ A with A = ±; 
this value is suggested by the result (4.42) obtained in the preceding subsection. 
The fits are very good, even for L m i n = 16 (see Table 14). Similarly good fits are 
obtained also with A = j, 1, 2. Finally, we also tried a logarithmic fit as in the 
4-state Potts model. The results (marked with L in Table 14) are also excellent. It 
is worth mentioning that the estimates of and x* T vary slightly from one type of 
fit to another, in some cases by several standard deviations. 

4.5 Dynamic quantities 

In this section we analyze both the integrated autocorrelation times T int ,o and 
the exponential autocorrelation times T exPj £>. Using standard power-law fits to the 
Ansatze T int ,o = AL Zint '° and T exPj £> = AI Zex P'°, we can extract the dynamic critical 
exponents. 21 However, the existence of multiplicative logarithmic corrections to the 
specific heat for the Ising and 4-state Potts models suggests, in view of the Li-Sokal 
bound, that similar multiplicative logarithmic corrections might occur also in the 
autocorrelation times. We shall look for such corrections in two ways: 

(a) by fitting to an explicit logarithmic Ansatz r = AL Z (log L)' p ; and 

(b) by studying the ratio r /Ch (for the X2 and ZF models we consider the ratio 

T I Cff,max)- 

4.5.1 Integrated autocorrelation time: power-law fits 

From Tables 4, 6 and 8, we see that the integrated autocorrelation time for 
Ai 2 (resp. Ai^ and Ai 2 T ) is always slightly smaller than the integrated autocorre- 
lation time for £ (resp. and E aT ). Furthermore, we see that the ratio of those 
autocorrelation times is a constant (i.e., independent of L) within statistical errors: 

I^ML(Y) = 0.978 ±0.011 (4.45a) 

21 We emphasize that in general z; nt) e> need not be equal to z exPt o • However, in SW-type algorithms 
it does appear that the autocorrelation function of the energy is very close to a pure exponential, so 
that Ti nt) £/T exP) £ approaches a constant (in fact, a constant very close to 1) as L — ► oo. So in these 
algorithms we do empirically seem to have 2; n t,£ = ^exp,f • 
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T ' mtMl (ZF) = 0.944 ±0.018; T ' mtM ^ (ZF) = 0.887 ± 0.020 (4.45b) 

^^(X2) = 0.857 ±0.018; T ' int > M ^ (X2) = 0.803 ± 0.017 (4.45c) 

Likewise, if we look at the ratio T- m t,£ aT I T^t^^ for the points ZF and X2, we see that 
in both cases that ratio is also consistent with a constant: 

Tint > £ -(ZF) = 0.960 ±0.017 (4.46a) 



^"int^o-T- 



(X2) = 0.962 ±0.013 (4.46b) 



It therefore suffices to consider the critical behavior of T int ,£ (for Potts) and Tjnt^ 
(for ZF and X2); all other quantities will have the same dynamic critical exponent. 

Potts model. If we fit the data from the 4-state Potts model point to a pure 
power-law function T int ,£ = AL Zint ' £ , we obtain a good fit already for L m i n = 16 (see 
Table 10). However, there seems to be a weak upward trend with L m i n} so to be 
conservative we choose L m i n = 32 as our preferred fit: 

2int,f(P) = 0.876 ±0.012 (4.47) 

with x 2 = 3.16 (4 DF, level = 53%). Notice that this value of z- mt ,£ is strictly greater 
than the effective exponent (a/v)^ = 0.768 obtained by a pure power-law fit [see 
(4.34)]. This implies that the Li-Sokal bound (1.1) is satisfied, but it is apparently 
not sharp. We can also compare the power-law estimates z- m t,£ and ajv for each 
Lmin separately (see Table 10), and again conclude that the Li-Sokal bound is 
always satisfied but that it is not sharp. Note, however, that this effective exponent 
(zint,£ — oi/v) e f[ ~ 0.11 is consistent with the true behavior of t^^/Ch being either 
a small positive power or a logarithm. 

If we compare our results for the 4-state Potts model with the embedding algo- 
rithm (see Table 4) to those quoted in [9] (which correspond to the direct algorithm 
of Section 3.1), we see that the ratio of the two autocorrelation times is more or less 
constant within statistical errors, and that there is no systematic trend as L grows. 
For these reasons, we conclude that they are proportional in the limit L — > oo: 22 

^-direct 

int,£ ^ = L516± o.035 . (4.48) 



_embedd 
'int,£ 



In particular, the direct and embedding algorithms belong to the same dynamic 
universality class (as was of course to be expected). It follows from (4.48) that our 



22 It is therefore hardly surprising that the dynamic critical exponent reported in [9] for a pure 
power-law fit, 2; n t,£ = 0.87 ± 0.02, is virtually identical to our value (4.47). If we fit our data for 
L < 256 as in [9], we obtain (for L min = 64) z mt)£ = 0.900 ± 0.025 with X 2 = 0.63 (1 DF, level = 
43%). This result is also compatible within errors with the one reported in [9]. 
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algorithm is 150% as effective as the standard SW algorithm for this model in terms 
of autocorrelation times. However, our algorithm needs roughly twice the CPU time 
for a complete sweep over the lattice (there are twice the number of points as in 
the Potts formulation) 23 ; so, in the end, the embedding algorithm is about 25% less 
efficient than the standard SW algorithm for the 4-state Potts model at criticality. 

Our estimate of z- int ,£ for a pure power-law fit, 0.876 ± 0.012, is very close to 
the result z- in t,£,ic = 0.92 ± 0.01 found by Wiseman and Domany [24] for the single- 
cluster version of the direct algorithm. 24 ' 25 Thus, the 2D Potts model with q = 4 
conforms to the behavior found previously [6] for the 2D Potts models with q = 2, 3, 
in which z- in t,£,ic ~ ^int,£,sw- On the other hand, this almost-equality seems not to 
occur for Ising models in dimension d > 3 [12,13]. 

ZF model. For the ZF point (Table 6), the estimates for ^mt,^ are quite stable, 
giving for L mm = 32 

2int,£ u (ZF) = 0.733 ±0.014 (4.49) 

with x 2 = 1-48 (3 DF, level = 68%). Once again, the Li-Sokal bound holds but it 
is not quite sharp, as here ajv = |. 

X2 model. At the point X2 (Table 8), the corrections to scaling for Tjnt^ seem 
to be significant, just as they are for CH,ma,x (see Table 13). The estimates for ^mt,^ 
tend to be systematically smaller as L m i n grows (although this effect is only on 
the borderline of statistical significance, due to the rather large error bars). Our 
preferred fit uses L m i n = 128: 

2int,£ u (X2) = 0.477 ± 0.028 (4.50) 

with x 2 = 0.39 (1 DF, level = 53%). Once again, the Li-Sokal bound holds but 
is not quite sharp, as a/v = 0.438 (our numerical estimate) or 0.4183... (exact 
value). 

Also here we can compare our result with the one found by Wiseman and Do- 
many [24] for the single-cluster algorithm. They obtain z- in t,£,ic = 0.61 ± 0.01. This 
value is very far from ours, even if we fit our data for 16 < L < 128 to facilitate the 
comparison with their data: in that case our preferred fit is z- mt ,£ = 0.553 ±0.012 for 
Lmin = 16 and % 2 = 0.32 (2 DF, level = 85%). Here the difference between z- in t,£,ic 
and z- int ,£,sw is nearly four standard deviations, so it seems that the two algorithms 
belong to different dynamic universality classes. 



23 We have directly measured the CPU ratio between these two algorithms and it is k, 1.9, i.e., 
very close to the guessed value of 2. 

24 We emphasize that 2; n t,f,ic is the dynamic critical exponent for the autocorrelation time of the 
single-cluster algorithm measured in units of "equivalent sweeps". This is d — j/u = 1/4 less than 
the dynamic critical exponent for the autocorrelation time measured in units of "cluster hits" . 

25 Actually, the Wiseman-Domany result was obtained with 16 < L < 128. If we fit our own data 
with this constraint, we get for L m i„ = 32 the value 2; n t,£ = 0.876±0.021 with \ 2 = 1-96 (1 DF, level 
= 16%). So our estimate of 2; n t,£ is not sensitive to L max ; and it differs from the Wiseman-Domany 
estimate of 2; n t,f ,ic by two standard deviations. 
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Ising model. Finally, we have reanalyzed the data of Baillie and Coddington [7] 
for the 2D Ising model (see Table 3). Our preferred fit is for L m i n = 100, giving 

2int,£(DIs) = 0.240 ± 0.004 (4.51) 

with x 2 = 1-67 (2 DF, level = 43%). In this case the Li-Sokal bound is clearly 
satisfied, as (a/v)^ 0.173 when L m i n = 128 (see Section 4.4.2); and once again 
the bound appears to be not quite sharp. 

In summary, the Li-Sokal bound is fulfilled in all cases, but apparently not as 
an equality: the effective exponents (z- in t,£ — «/^)eff range from 0.04 to 0.11. 
However, we know that in two of the four cases — namely, the Ising model and the 
4-state Potts model — the leading behavior of the specific heat is not merely a power 
law, but rather contains multiplicative logarithmic corrections. So the question of 
the sharpness of the bound cannot be answered until we take into account the exact 
leading term of the specific heat. This will be done in the next subsection. 

4.5.2 Might the Li— Sokal bound be sharp modulo a logarithm? 

It is well known that the true leading behavior of the specific heat is not given 
by a/i/ = 0.768 (or any other power law) for the 4-state Potts model, nor by 
a/v = 0.173 (or any other power law) for the Ising model. Rather, the behavior is 
-Llog~ 3 / 2 L and logL, respectively. A similar problem arises for the autocorrelation 
times: to answer the question of the sharpness of the Li-Sokal bound, it is necessary 
to guess the "exact" leading behavior of T int ,£- In analogy with Ch , we may entertain 
the possibility that T int ,£ contains a multiplicative logarithmic term (at least for the 
Ising and 4-state Potts models). 

Our first approach is to consider the ratio T int ,£ /Ch- 26 For all four models, we 
find that this ratio is an increasing function of the lattice size L. We then distinguish 
three possible asymptotic behaviors: If t^^/Ch tends to a constant as L — > oo, 
then the Li-Sokal bound (1.1) is sharp; if it grows like some logarithmic function, 
then the bound is sharp modulo a logarithm; and if it grows like some power law, 
then the bound is not sharp. It is reasonable to hope that the ratio r mtj £ /C H will be 
less affected by corrections to scaling than either Ch or T int ,£ separately. (Of course, 
this hope may also be false!) For the Ising model we took T int ,£ from Ref. [7] and 
the specific heat from the exact finite- volume solution given in Ref. [53]. For the 
other three models, we used our numerical data; in computing the error bar on the 
ratio, we used the triangle inequality, which of course yields only an upper bound 
on the true error bar. This overestimate of the error bars in the three non-Ising 
cases should be taken into account when interpreting the results. 

We tried to fit the ratio T int ,£/ 'Ch to various different Ansatze (see Table 16). The 
first is a pure power-law function. In all cases the fit is very good and the estimates 
of the power are very small (between 0.05 and 0.12); the power seems to increase 

26 More precisely, we consider Ti^s/Ch for the Ising and 4-state Potts models, and Ti n t t £^/CH,max 
for the X2 and ZF models. 
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slightly as we go from the Ising model to the 4-state Potts model. Note also that the 
fit for the Ising model requires L m i n = 100 in order to get a reasonable % 2 , while for 
the other three models an excellent % 2 is obtained already for L m { n = 16. This arises 
from the very accurate data in the Ising case, which permit the observation of very 
small corrections to scaling, in contrast to the rather larger (and overestimated!) 
error bars in the three non-Ising cases. 

Often such small powers indicate that the true behavior is logarithmic. Indeed, 
a logarithm can be very well mimicked by a power law over any not-too-wide range 
of L. We therefore tried various combinations of logarithmic functions. The first 
one was a function A\og p L. The quality of fit for the Ising model is rather inferior 
to that obtained with a power-law function, although the confidence level is still 
reasonable (27%). However, for the three other models, the fits are very stable and 
give excellent values of the % 2 . Next we tried a fit to A + BlogL; this gives very 
good results for all four models. 

Finally, we have tried to check whether the ratio T int ,£/ 'Ch approaches a constant 
as L — > oo. First we tried a fit to a pure constant A. Even by eye one can see 
that the values of t^^/Ch for different L exhibit statistically significant deviations 
(e.g., consecutive values differing by at least two standard deviations) for L < 
128, even for the non-Ising models where we already know that the error bars are 
overestimated. So, not surprisingly, it is impossible to get a decent fit to a constant 
with L m i n < 128. However, for the non-Ising models with L m i n > 256, the ratios 
7"int,£ /Ch are consistent with a constant A, at least if one takes at face value the 
overestimated error bars: see Table 16. It follows that no reliable information can 
be obtained on the manner of convergence to a constant for these models, if one 
takes L m i n > 256. Next we tried fitting the ratio t-^^/Ch to A + BL~ A with 
A = 2, 1, |, |, and |, and also to A + B/logL ("A = X log"). For all the models, 
the fits with A = 1,2 are both implausible and unreliable, as the parameter B 
becomes very large (\B\ > 10) and does not stabilize as L m i n grows. For the Ising 
model, only the fits with A = |, | have reasonable % 2 values (confidence levels of 
14% and 28%, respectively). Note, in particular, that the logarithmic-correction 
Ansatz gives a poor fit (level = 5%). For the 4-state Potts model we always obtain 
unreasonably large values of \B\ (> 10), so we cannot trust these fits, even if they 
have reasonable values of % 2 . We therefore rule out this scenario for the 4-state 
Potts model. The good values of the % 2 are surely due to the overestimated error 
bars. This is a warning against trusting the % 2 values for the other two non-Ising 
models. Finally, we get reasonable fits for the X2 and ZF models for < A < |; 
but perhaps the good % 2 values are not to be trusted. 

Let us discuss the above fits. For the Ising model we have been able to associate a 
reliable error bar to the ratio t-^^/Ch-, so we can trust the % 2 values. We conclude 
that an asymptotically bounded ratio with additive corrections like either A = 
X log or A > | is very unlikely to occur. The most plausible scenarios (i.e., 
highest confidence levels) are the pure power law with p = 0.060 ± 0.004 or a simple 
logarithmic growth A + BlogL. However, a logarithmic power-law behavior with 
p = 0.315 ±0.020, or an asymptotically bounded behavior with additive corrections 
given by | < A < |, cannot be completely ruled out. 
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For the other three models we cannot trust the % 2 values, as the error bars are 
overestimated. We have followed three criteria to interpret our results and decide 
which are the "best" fits: 

• Absolute \ 2 value. If the \ 2 value of a given fit is not good, then the fit is 
surely poor, as all the error bars are overestimated. 

• Relative x 2 values. If two fits to the same data exhibit vastly different values 
of % 2 /DF, then the fit with the larger % 2 /DF can be considered less plausible 
(all other things being equal). 

• Reasonable results. When fitting to a constant plus additive corrections (A + 
BL~ A ), we expect that the parameter B should remain not too large (e.g., 
\B\ ^ 10). If this does not occur, we tend to distrust the fit. 

• Value of L m i n . If L m i n is taken large enough (i.e., > 256), then we can fit 
the data for the three non-Ising models with almost any additive-correction 
Ansatz, as all the values of t^^/Ch are already consistent with a constant 
within the (overestimated) errors. Thus, given two fits with similar % 2 /DF 
values, we tend to trust more the one with smaller L m i n . 

The 4-state Potts model is the clearest case. The asymptotic constancy of 
7"int,£ /Ch is clearly ruled out due to unreasonably large parameter B in all the fits. 
Moreover, all the asymptotic-constant fits need L m i n = 64 to obtain a decent % 2 
(even with the overestimated error bars) whereas L m i n = 16 is sufficient for the pure 
power-law and A + B log L scenarios. Even the logarithmic power-law Ansatz needs 
Lmin = 32 to achieve a comparably good % 2 . If we consider the % 2 when L m i n = 16, 
we see that the best fits have % 2 = 1.39 (power-law) and 1.98 (A + BlogL) } while 
the others have 3.47 (A = §), 3.64 (Alog p L) and 5.61 (A = ±). We conclude that 
the two most likely scenarios are the pure power-law with p = 0.018 ± 0.012 and 
the simple logarithmic behavior A + BlogL. 

For the X2 model we arrive at the same conclusion. The best fits need L m i n = 16, 
and the rest need at least L m i n = 32. If we consider the % 2 when L m i n = 16, we 
see that the best fits have % 2 = 1.28 (power-law) and 1.43 (A + .BlogL), while the 
others have \ 2 = 1.96 (A = |), 2.39 (logarithmic power-law) and 2.68 (A = |). 

Finally, the ZF model is the least clear case. Here most of the fits are reasonable 
with L m i n = 16. Thus, we cannot decide among the pure-power scenario, the two 
logarithmic scenarios, and the scenario of asymptotic constancy with corrections to 
scaling given by | < A < \. 

On theoretical grounds one might expect some kind of "continuity" in the be- 
havior of the ratio t^^/Ch along the self-dual curve: that is, one might expect the 
same scenario to hold everywhere along the curve (except perhaps at the Ising and 
4-state-Potts points, where there might be additional logarithmic effects). There 
are only two scenarios that are consistent with the data in all four cases: 





(4.52) 
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Using the theoretically known exact behavior of C#, we can investigate the validity 
of the Ansatze (4.52) directly on T int ,£- For the Ising model these Ansatze become 



A'L P log L 
A' log L + B' log 2 L 



r int , £ (DIs) « { 4/i r 6 | n/i 2r (4.53) 



A fit to the first Ansatz with L mm = 100 has x 2 = 1-35 (2 DF, level = 51%) and 
gives p = 0.051 ± 0.004, while a fit to the second Ansatz with L m i n = 100 has 
X 2 = 1.69 (2 DF, level = 43%). For the X2 model the Ansatze are 



(A + B log L)L° 



r mt , f .(X2) « \ /a _ „ 1 ... rNr0 . 4183 (4.54) 



The first (pure power-law) Ansatz has already been studied in the preceding sub- 
section, yielding x 2 = 0.39 (1 DF, level = 53%) with L mm = 128 [see (4.50)]; the 
second Ansatz yields x 2 = 0.42 (1 DF, level = 52%), also with L m i n = 128. So 
both fits are good, and there is nothing to distinguish them. For the ZF model, the 
Ansatze are 

f A l L P+2/3 

W,(ZF) « ( (A + Blog£)ri/ 3 <« 5 > 

The first (pure power-law) Ansatz has already been studied in the preceding sub- 
section, yielding x 2 = 1-48 (3 DF, level = 68%) with L mm = 32 [see (4.49)]; the 
second Ansatz yields x 2 = 1-53 (4 DF, level = 82%) with L m i n = 16. So in this 
case there is a slight preference for the logarithmic Ansatz. Finally, for the 4-state 
Potts model the Ansatze are 

, m [ A'L P+1 log~ 3 ^ 2 L /Arr , 

(P) * Ua + b^dl^l (4 ' 56) 

The first one gives the power p = 0.13 ± 0.03 for L m i n = 128 with x 2 = 0.96 (2 DF, 
level = 62%). The second one yields x 2 = 0.98 for the same L m i n (2 DF, level = 
61%). Finally, one can try the fit 

r int , f (P) « ALlog- p ' L, (4.57) 

in which one imposes Zi nt ,£ = 1 X log~ p and attempts to find the multiplicative 
logarithmic exponent p' . Here the stability of the results is not very good (see 
Table 17). One could argue that the fit with L m i n = 16 already has a reasonable 
X 2 value, but we are inclined to be conservative and prefer the fit with L m i n = 128 
(which of course has much larger error bars): 

p'(P) = 0.776 ± 0.180 (4.58) 

with x 2 = 1-09 (2 DF, level = 58%). This result is compatible (within 2 standard 
deviations) with the expected value of p' = |. If we take into account the value 
reported in Table 16 for the fit T- m t,e/CH = A\og p L [p = 0.543 ± 0.073], we see 
that our estimates for p and p 1 are compatible within errors (i.e., p' = | — p = 
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0.957 ± 0.073). However, that value of p is not compatible with the value of p' = | 
corresponding to Ansatz (4.56b). Actually, our estimate of p' lies midway between 
the Ansatz and the value coming from our estimate of p. 

From the above results it is difficult to tell which is the true asymptotic behavior 
of Tint.f/Cff (and thus of Ti ntj ,f). To disentangle this we would need much larger 
lattices, as well as much better statistics on the lattices L > 128. 

Remark. Let us also test the conjecture proposed by Baillie and Coddington 
[7] for the Ising model: 

T mti£ ^ (A + B log L)L^ V , (4.59) 

where f3 is the static critical exponent for the spontaneous magnetization; note that 
/3/v = | everywhere on the AT self-dual curve. For the Ising model, the fit is fairly 
good for L m i n = 100, giving x 2 = 0.50 (2 DF, level = 78%). The same conclusion 
applies to the X2 model (L mm = 128, x 2 = 0.42, 1 DF, level = 52%) and to the ZF 
model (L m i n = 128, x 2 = 0.44, 1 DF, level = 51%). However, the fit is rather poor 
for the 4-state Potts model: for L m i n = 256 we only get x 2 = 4.43 (1 DF, level = 
4%). 

4.5.3 Exponential autocorrelation time 

The exponential autocorrelation time is for an observable O is defined as 27 

T exp ,o = hm — ■ : — . (4.60) 

t^oo - log \poo(t)\ 

This autocorrelation time measures the decay rate of the "slowest mode" of the 
system, provided that this mode is not orthogonal to O. 

The critical behavior of T exPj £> is, in general, different from the behavior of T int ,o- 
This fact can be seen from the standard dynamic finite-size-scaling Ansatz for the 
autocorrelation function pooif)'- 

Poo(t; L) « \t\-*°ho ( — ; ^) • (4.61) 

(Here the dependence on the coupling constants has been suppressed for notational 
simplicity.) Summing (4.61) over t, it follows that 

Tint,0 ~ , (4.62) 

or equivalently, 

Zint,0 = (1 - P0)Zexp,0 ■ (4.63) 

Thus, only when po = do we have z- in t t o = Zex P ,o [1 5 4]. 

Here we consider the exponential autocorrelation time of the energy £ for the 
4-state Potts model and of the observable for the X2 and ZF models. We expect 



27 Strictly speaking, the "lim" should be replaced by "limsup", as in (4.4). But in virtually all 
practical applications, the limit really exists. 
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that a similar behavior would be found for £ ar for the X2 and ZF models, and for 
the squared magnetizations in all three models; but it did not seem worthwhile to 
us to carry out this analysis in detail. 

If we plot the estimated autocorrelation function pssit), we see that it fits beau- 
tifully to a pure exponential for t > T- m t^/2: as examples, see Figures 2 and 3, 
which represent the data for the 4-state Potts model with L = 16, 32 and L = 256, 
respectively. So it makes sense to extract estimates of r exPj £ from the tail of the 
autocorrelation function. More precisely, for each run we estimated r exPj £ by fitting 
log p ££ (t) = -A - Bt for the interval t mm <t< t max ; obviously r exPi £ = l/B. By 
studying the goodness of fit (i.e., the % 2 value and the corresponding confidence 
level) as a function of i m ,- n and t maX} we can choose the "best" fit. 

This fit is, however, extremely subtle, because the Monte Carlo estimates of 
Pee{t) for different t are in general (highly) correlated. The full covariance matrix 
C p for these random variables can in principle be computed using the autocorrelation 
function itself, at least in the approximation that neglects the fourth cumulant of the 
stochastic process: see equation (B.l) in Appendix B. If one assumes for simplicity 
that the autocorrelation function is a pure exponential, then the formula simplifies 
further: see equation (B.2). From this formula we already see that the off-diagonal 
terms in C p are comparable in magnitude to the diagonal ones; therefore, it is 
unlikely to be sensible to neglect them. 

Nevertheless, one may try the crude approximation in which all off-diagonal 
terms in C p are dropped, and see what happens. Using the self-consistent procedure 
explained in detail in Appendix B, we obtained both an estimate of r exPj £ (along with 
its error bar) and the error bars for the autocorrelation function pee{t)- For almost 
all values of i m ,- n and t maX} we obtained a perfect fit (confidence level 100%), with 
values of the % 2 much smaller than the number of degrees of freedom. However, 
the variation of the estimates of r exPj £ as a function of i m ,- n and t max was much 
larger than the error bars given by the fit. This indicates that the error bar for the 
estimated r exPj £ was not correctly computed. Clearly, the neglect of the off-diagonal 
covariances is unjustified, as we already expected on theoretical grounds. 

We therefore redid the fit using the full covariance matrix C p} again using the 
self-consistent procedure described in Appendix B. As before, we systematically 
varied i m ,- n and t max . The dependence of the results on t max is usually slight; but 
the dependence on i m ,- n is moderately strong. We report our results in Tables 18-19 
for the 4-state Potts model 28 , in Tables 20-21 for the X2 model, and in Tables 22-23 
for the ZF model. For each model we have presented two different tables, one with 
tmax = 4r intj£ and the other with t max = 3r intj£ . 

In all cases we have a bad fit (level <C 1%) whenever i m ,- n <C 0.5ri ntj £ (this 
happens irrespective of the value of t max ). As an example of such a bad fit, we show 
in each table the fit with i m ,- n = 1. Clearly, the energy autocorrelation function is 
significantly different from a pure exponential for very small t. However, we obtain 
reasonable % 2 values as soon as i m ,- n > 0.5ri ntj £, indicating that the autocorrelation 
function becomes very close to a pure exponential for t > 0.5ri ntj £. 

28 The results for L = 1024 have not been quoted, as the statistic is rather poor (~ 1500ri nt) £). 
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For the 4-state Potts model (Tables 18-19) we have rather stable results for 
L = 16, 32, 64. In most of these fits the confidence level remains reasonable, so 
we think we can fully trust these results. However, for L > 128 we begin to have 
difficulties: the self-consistent procedure does not always converge, and the results 
with t max = 4:Ti nt ,£ sometimes differ from those with t max = 3ri ntj £ by about two 
standard deviations. Furthermore, for these large lattices we consistently obtain 
unusually high confidence levels (mostly > 99%); we do not understand why this 
happens. For these lattice sizes the troubles seem to be somewhat less severe when 
tmax = 3ri ntj £. So we tend to trust better these latter fits when L > 128. 

For the X2 model (Tables 20-21) both the stability of the results and the good- 
ness of fit are excellent for L = 16 and L = 32. For L = 64 we have bad fits when 
tmax = 4Ti nt) £ a , (level < 1%), and somewhat better ones when t max = 3^^^^ (level 
~ 8%). Also for L = 128 we obtain somewhat more stable and consistent results 
when t max = Srint^. For L = 256 the fits are rather poorer; for L = 512 they are 
good, but the results for the two different values of t max differ by three standard 
deviations. 

Finally, for the ZF model (Tables 22-23) we have good stability and goodness 
of fit for L = 16 and L = 32. For L = 64 the fits are usually very poor (level < 
5-8%), but they are generally less bad (and also more stable) for t max = 37"^^ 
than for t max = 47i n t ) £ a ,. For L > 128 we also find the results with t max = 37"^^ 
somewhat better than for t max = 47"^^ ; again we find inexplicably high confidence 
levels (often > 99%). 

From the above discussion we can already see that we have good and stable fits 
only for the smaller lattices (L < 32 and in some cases L = 64), where the statistics 
are better. Apparently, in order to have a decent estimate of r exPj £ we need a run 
length of at least 6 X 10 4 T int ,£ and possibly more. 

To decide which are the "best" fits we use the following criteria: we choose the 
largest interval [t m i n} t max ] such that 

(a) the x 2 value is reasonable (e.g., confidence level > 10%); and 

(b) there is some consistency (within error bars) with the values coming from 
"nearby" choices of i m ,- n and t max . 

In Table 24 we present what we consider to be the "best" fits for each model and 
each lattice size L. We include the interval [t m i n} t max ] and the ratio ri ntj £/r exPj £. 
The error bar on this ratio was computed using the triangle inequality, as we do 
not know what is the covariance between our estimates of T int ,£ and r exPj £. We think 
these estimates are reasonably reliable for L = 16 and 32, somewhat reliable for 
L = 64, and possibly unreliable for L > 128. 

From Table 24 we see that for each model, the values of the ratio T int ,£ / T eX p,£ are 
more or less constant when L < 64. If we fit the ratios for 16 < L < 64 to a constant, 
we obtain reasonable fits for the three models: T int ,£ / T eX p,£ = 0.936 ± 0.015 for the 
4-state Potts model ( X 2 = 1.50, 2 DF, level = 47%), r intj ^/r exPi ^ = 0.966 ± 0.015 
for the ZF model ( X 2 = 0.10, 2 DF, level = 95%) and r int ',^/r exPi ^ = 0.980 ± 0.013 
for the X2 model (x 2 = 0.50, 2 DF, level = 78%). However, for the 4-state Potts 
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and ZF models, this ratio seems to decrease when L > 128. Now, these points are 
precisely those where we had troubles obtaining the value of r exPj £, so this decrease 
might be due simply to a poor determination of the exponential autocorrelation 
time; and it is in any case only about two standard deviations. 29 On the other 
hand, this decrease might indicate that T int ,£ / T eX p,£ is tending to zero as L — > oo; 
in this case the dynamic critical exponents z- mt ,£ and z exPj £ would be different (the 
latter would be larger), and the exponent p£ in (4.61) would be strictly positive. 

The result for the X2 model, by contrast, is completely consistent with a constant 
ratio Tjnt^ /rexp,^ , indicating that ps w = 0. This fit to a constant, using all the data 
(16 < L < 512),' gives 

^^(X2) = 0.976 ±0.011 (4.64) 

with x 2 = 2.14 (5 DF, level = 83%). 

Due to the ambiguities in the determination of r exPj £ for L > 128 in all the 
models, we are unable to come to any definitive conclusion on whether z- mt ,£ = Zex P ,£- 
But it does seem likely. 



5 Conclusions 

We have performed a high-precision Monte Carlo study of the symmetric Ashkin- 
Teller model at several points on the self-dual (critical) curve, using a Swendsen- 
Wang-type algorithm. We have considered both the static behavior of the models 
(known exactly) and the dynamic behavior of the algorithm. 

We have had great difficulties in obtaining the correct leading behavior whenever 
this is not simply a power law plus additive power-law corrections. These difficulties 
occurred both for static quantities (specific heat in the Ising and 4-state Potts 
models) and for dynamic quantities (autocorrelation times in all models). Unless 
we have some theoretical input, it is almost impossible to distinguish between power- 
law and logarithmic behaviors when the range of lattice sizes L is not extremely 
large (in our case, 16 < L < 1024). 

This issue makes it very problematic to tell whether the Li-Sokal bound (1.1)/ (1.2) 
is sharp or not. Our results seem to indicate that there are only two likely scenarios: 
the Li-Sokal bound fails to be sharp either by a small power [i.e., t^^/Ch ~ AL P 
with 0.05 ^ p ^ 0.12] or by only a logarithm [e.g., t^^/Ch ~ A- + -BlogL]. Either 
one of these scenarios is consistent with our data at all four points on the AT self- 
dual curve. Larger lattices and much better statistics will be needed to distinguish 
between them. 

We have also presented a new method for estimating the exponential autocor- 
relation time, which takes into account the full covariance matrix for the sample 

29 If we fit all the data (16 < L < 512) we get the following ratios: Ti Dt t £ / T exp £ = 0.914 ± 0.013 
for the 4-state Potts model ( X 2 = 13.84, 5 DF, level = 2%) and r mt) ^/r e ' xP) ^ = 0.945 ± 0.014 for 
the ZF model (\ 2 = 10.84, 5 DF, level = 5%). These confidence levels on the order of 5% are as 
expected from the two-standard-deviation discrepancies. 
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autocorrelation function. To do so is essential to obtain reliable results, as the val- 
ues of the sample autocorrelation function are strongly positively correlated. The 
quality of the estimates of T exPj e> depends strongly on the accuracy of the available 
data: we seem to get reliable estimates of T exPj e> only when the run length is at least 
« 60000r intiO . 

A Proof of Li— Sokal bound for the direct Ashkin— 
Teller algorithm 

We can easily extend the proof of the Li-Sokal bound for the g-state Potts model 
[9] to the direct algorithm for the AT model (defined in Section 3.1). Also in this 
latter case, the transition matrix can be written as a product 

PSW = -Pbond-Pspin , (A.l) 

where Pbond (the update of the bond variables) and P sp in (the update of the spin vari- 
ables) are given by the conditional expectation operators E( ■ |{cr, r}) and E( ■ |{m, n}), 
respectively. 

As in [9], we are going to compute explicitly the autocorrelation function at 
time lags and 1 for several bond densities. Then, using some general properties of 
reversible Markov chains, we will deduce lower bounds for the autocorrelation times 
T- m t,A (for certain observables A) and r exp . These will in turn imply lower bounds 
on the dynamic critical exponents z- mt ,A and z exp . 

Let us consider the bond occupations 30 

(A. 2a) 

(xy) 

M = J2 n -y ( A - 2b ) 

(xy) 

O = £ m xy n xy (A. 2c) 

(xy) 

We will follow the notation of [9] and henceforth write the Kronecker deltas for a 
bond b = (xy) as 8 ab = 8 ax)(Ty and S n = S TxjTy . 

From (3.6) we can read off the expectation value of the bond variable m& con- 
ditional on the spin configuration {cr, r}: it is 

E{m b \{a,T}) = q 1 8 ab = {p 1 + p 2 ) 6 ab , (A. 3) 

where qi } pi and p 2 are the probabilities appearing in Steps la and lb of the direct 
algorithm (Section 3.1). 31 The important fact here is that qi = pi + p 2 : this means 



30 Do not confuse this M with a magnetization! 

31 In this appendix we are assuming that the system is homogeneous, i.e., that the couplings J, 
J' and K do not vary from bond to bond. An inhomogeneous system can be treated by an obvious 
generalization. 
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that E(m b \{a } r}) does not depend on the r configuration. Likewise we have 



E(n b \{a,T}) = r 1 S Tb = (pi+p 3 )S Tb 
E(m b n b \{a, r}) = p l 6 !Tb 6 Tb 

The other conditional probabilities we need are 



E(m b m b i 
E(n b n b , 
E(m b n b i 
E(m b n b m b i 
E(m b n b n b i 
E(m b n b m b in b i 



q 2 S ab S ab , iib^b' 

q t 8 ab if b = b' 

r\6 Tb 6 Tbl iib^V 
r\ 6 



if b 



qin S ab S Tb , lib^V 
pi S ab S Tb if b = V 



I Pi S ab S Tb if b = V 



f Pl^l &a b &T b &T h , if 6 7^ 6' 



I Pi ^(7 b ^ 



6^6 



if 6 = 6' 



Pi $a b & 



b^T b 



if 6 



(A.4a) 
(A.4b) 



(A. 5a) 
(A. 5b) 
(A. 5c) 
(A.5d) 
(A.5e) 
(A.5f) 



These simply reflect the fact that in Step 1 each bond is updated independently of 
each other bond. 

From (A. 3) and (A. 5a) it is easy to compute the mean values (Ai) and (Ai 2 ) } 
and hence vax(M) = (M 2 ) - (M) 2 : 



(M) 
(M 2 ) 
vax(M) 



2 y|(l+K) 

?1 2 /C2\ 



4 
2V 



(£ 2 ) + qtV'il + 2E a ) + qi (l - qi )V(l + E a ) 



— C H ,<7<7 H ^ (1 + K) 



(A. 6a) 
(A. 6b) 
(A. 6c) 



where E a = (1/2V)(£ (J ) is one of the energies defined in Section 4.2, and Cn,aa = 
(l/2V)var(£ CT ) is the corresponding specific heat. The unnormalized autocorrelation 
function of Ai at time lag is precisely Cmm{§) = var(A^). 

The corresponding autocorrelation function at time lag 1 is given by 

Cmm(I) = (M(O)M(l)) - (M) 2 = ™r(E(M\{vT})) = var(P bond X) . (A.7) 

Now, the energy-like operator PbondA'f is equal to 



P hond M = E(M\{*,t}) 



This implies that 



-var 



(&) 



-(2V + £ a ) 



2V^C H , aa . 



(A.8) 



(A.9) 
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Thus, the normalized autocorrelation function for the operator Ai at time lag 1 is 
given by 



n m _ C MM (1) _ 2(l- gi )(l+iS g ) 

} " C^(0) - qi C H ^ + 2(1 - 50(1 + E a ) ■ (A - 10) 

Note now that when we approach any point on the critical curve, the quantity 
qi = 1 — e~ 2 ( J+A ) remains positive and less than 1 [i.e., q\ — s- qi )C -rit £ (0, 1)], while 
the energy E a remains greater than —1 [i.e., E a — > E afir a > — 1]. It follows that 

, . const 
Pmm(1) > I-7; (A.ll) 

uniformly in a neighborhood of that critical point. 

The rest of the argument given in [9] can now be transcribed verbatim. The cor- 
relation functions of Ai under Psw are the same as under the positive-semidefinite 
self-adjoint operator P' sw = -Pspin-Pbond-Pspin- This fact implies that we have a spec- 
tral representation 

PMM(t)= I' \ ltl du(\) (A.12) 
Jo 

with a positive measure dv. From this equation we conclude that 

PMM(t) > PMM(lf l ■ (A.13) 

Using the definition (4.3) of the integrated autocorrelation time and the definition 
(4.4) of the exponential autocorrelation time, we arrive at the following bounds 

1 1 + pmm(1) ^ , ri ( k ~\ a \ 

Tint,X > - 7TT > COnst X C H ,aa (A. 14a) 

2 1 - Pmm(i-) 

T exp ,M > i — TTT > Const X C H ,oo (A. 14b) 

log Pmm\L) 

Let us now assume that the autocorrelation times diverge (for a finite system of 
size L at criticality) as L Zint ' M and L Zex p< M , respectively, and that the matrix element 
Ch,uu of the specific-heat matrix diverges as L a l v . We then conclude that 

OL 

Zint,M^ Z exp,M > — , (A. 15) 

which is the result of Ref . [9] . 

The only way that the bound (A. 15) could fail is in case the matrix element Cn,aa 
fails to diverge as L a l v (e.g., by diverging with a smaller power or by being bounded). 
Now, the matrix Ch does have an exactly marginal eigenvalue everywhere on the 
self-dual curve (2.8) of the symmetric AT model, so that the component of Ch 
tangent to this self-dual curve (namely, CH,mm) is bounded as L — > oo. 32 However, 



32 More generally, in the non-symmetric AT model on the self-dual manifold (2.7), there are two 
marginal eigenvalues, corresponding to the two tangent vectors to (2.7). 
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this marginal direction is never exactly cr, so the preceding proof does in fact always 
succeed: Cn,aa does always diverge as L a ^ . Nevertheless, some extra generality — 
as well as slightly sharper constants in the lower bound (A. 14b) on r exp — can be 
obtained by choosing in place of Ai a more general bond observable 



X = Cl M+c 2 N + c 3 (D 



(A.16) 



where ci, c 2 and c 3 are arbitrary real constants. This case trivially includes the 
preceding one. Using the techniques described above and after some algebra we 
arrive at the following results 



(X) 

C xx (0) 
C xx (l) 



2V 



ClTl (1 + E T ) + ^-(1 + E a ) + ^T"(l + E T ) 



2 v 2 

2V (C H ,p bond xp boIld x + E x ) 
2V Cff : p bond xp bond x 
E x 



1 



^ H,P ho - ad XP ho - ad X ~\~ E x 



where we have defined the following quantities 
Et = E T + E a + E ar 
P hond X = E(X\{a,r}) 



(A. 17a) 

(A.17b) 
(A.17c) 

(A.17d) 



(A. 18a) 



+ 



c 3 Pi 



+2V 



4 J " ' V 2 4 / 
2ciri + 2c 2 qi + c 3 pi 



C "<TT 



(A. 18b) 



^-ff i -Pb o n d X Pb o n d X 



E 



x 



—v&r(P hond X) 
1 + Et 

[c|pi(l - Pi) + 2c 1 c 2 (p 1 - q\) 



(A. 18c) 



+ 



+ 2cic 3 pi(l - ri) + 2c 2 c 3 pi(l - 9i )] , (A.18d) 



and the rest of the observable quantities are defined in Section 4. 

Using the parameters (ci, c 2 , c 3 ) we can choose the energy-like operator Pbond^ hi 
such a way that it contains a non-zero projection on the most divergent eigenvector 
of the matrix Ch- This is always possible as we have enough freedom. On the other 
hand, the term E x should be bounded from below by a strictly positive number. 

Let us analyze an interesting particular case: the symmetric model (r^ = ^) on 
the self-dual curve. The natural choice is c\ = c 2 and the previous formulae can be 
simplified to 



-Pbond^ 



+ Cl?l ) Ouj + —^-tar + 2V ( Ci^i + — ^— 



(A.19a) 
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E x = c{ qi (l - qi )(l + E u ) 

+ ^%|pi(l - Pi) + 4c lC3 pi(l - ?i) + 2cl( Pl - ql)] (A.19b) 

We can choose any ratio c 3 /ci such that Pbond^ is not a multiple of <-t m i n = £ LU -\-a£ (JT} 
where a is the parameter defined in (4.26). This means that on the self-dual curve 
the ratio c 3 /ci could be anything different from — |sinh2J. In particular, we can 
make the choice (ci,c 3 ) = (1,0), leading to Pbond^ = 2<1i£lu + const and to the 
bound 

Pxx(X) = 1 ~ , r E \ p (A-20a) 

2 

E x = qi (l- qi )(l + E u ) + ^^(l + E T ) (A.20b) 

Let us now restrict attention to the part of the self-dual curve for which the specific 
heat is divergent, namely the interval between the decoupled Ising point (DIs) and 
the 4-state Potts point (P). In this interval we have J } K > 0, so Griffiths' first 
inequality implies that the energies E w and E ar are > (and hence so is Ej). 
Moreover, both qi and pi — q\ are strictly positive on this interval. Since the 
specific heat Ch,ww is divergent everywhere on this interval, the proof of the bound 
(1.2) is complete. 



B Fitting (highly correlated) autocorrelation func- 
tions 



Let p(t) be the normalized sample autocorrelation function for some particular 
observable, measured in a Monte Carlo run of length N. These measured values 
t=~-(N -l) are °f cours e random variables; as such they have a covariance matrix 
C p , which is given by a standard formula from time-series analysis [51]: 



cov 



[p(r), p(s)} 



+ 00 



1 E 

N „^ 



p(m)p(m + r — s) + p(m + s)p(m — r) + 2p(r)p(s)p(my 



2p(r)p(m)p(m — s) — 2p(s)p(m)p(m — r)j + O 



N 2 



(B.l) 



where {p(t)} are the true values of the autocorrelation function. This formula is 
valid in the limit N — > oo, provided that the stochastic process is Gaussian. [If 
the process is not Gaussian, then we have to include also terms proportional to the 
fourth cumulant /c 4 (m, r, r — s).] In our case the stochastic process is of course not 
Gaussian, but we may hope that it is "not too far from Gaussian". So let us for 
simplicity take formula (B.l) as correct. 

The qualitative import of (B.l) can be understood by examining the special case 



of a pure exponential decay p(t) = e 



-\t\h 



In this case (B.l) reduces to 



cov 



[p(r),p{s)} 



1 

~N 



e -\rs\lr , 



A + 



1 + e 



-2/r' 



-2/r 
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-(r + s)/r ( r + s + 



1 - fi-2/r 



+ 0[±) .(B.2) 



We thus see that the off-diagonal terms in C p (i.e., r ^ s) are comparable in 
magnitude to the diagonal ones (r = s). In other words, the sample autocorrelations 
p(t) for different time lags t are strongly positively correlated , and any valid analysis 
method must take proper account of this correlation. 33 

Our Ansatz for the autocorrelation function p(t) will be the following: 

p(t) = \ m fOT °^<^ (B .3a) 

I Ae-'/^p for t > t mm 

p(t) = p(-t) (B.3b) 

where p(t) is the autocorrelation function at time lag t measured in the Monte Carlo 
simulation, and r exp will be chosen by least-squares fitting (see below); here i m ,- n 
is some chosen cut point, to take account of the fact that the behavior of p(t) for 
small t need not be exponential. Now, for each i m ,- n we can compute explicitly the 
sums appearing in (B.l), when p(t) is given by (B.3). Indeed, all the terms in (B.l) 
can be written in terms of p(t) and its convolution 

+00 

p(s) = ^ p(m)p(m — s) . (B-4) 

m = — 00 

The sum (B.4) can be split in two pieces: one piece contains only p(t) } s with t > t m8n , 
and the other piece contains the rest. The first piece can be summed exactly, giving 
the result 

s + fmm - l s + fmm - l g — (2t rnln + s) /r ex p 

K 5 ) = p(m)p(m-s)+ Z p(m)p(m - s) + 2A 2 _ 2/ ^ (B.5) 

m = l m = .s 1 e CXP 

Thus, given {p(t)} } t m i n} A and r exp we can compute the covariance matrix C p given 
by (B.l). With this matrix, we can perform the standard weighted least-squares fit 
[56] to log p(t) = a + bt for the interval i m ,- n < t < t max and obtain new estimates 
for A = e a and r exp = —1/6. We iterate this process until we reach a fixed point, 
for which the input and output values of A and r exp are equal. In practice, we 

initialized this self-consistent process by supposing that p(t) = e -l*l/ T -p with 

here the value of T int is of course our estimate from the Monte Carlo simulation, 
using the usual self-consistent truncation window of width 6ri nt [52, Appendix C]. 
We have followed this procedure both for the fits with the full covariance matrix 
and for those with only a diagonal covariance matrix. 



33 A similar problem arises in fitting the spatial correlation function to an asymptotic exponential 
decay in order to extract the mass gap [55]. 
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This procedure was implemented using MATHEMATICA, which allowed us to 
control accurately the numerical precision of the calculation. This is especially 
important when inverting the full covariance matrix, as in some cases we obtained 
nearly-singular matrices. In practice, this method converges well for the smaller 
lattices (the number of iterations needed is usually < 10). However, when the 
data are sufficiently poor (run length < 60000ri nt ), we noticed some cases of cyclic 
behavior instead of convergence to a single fixed point. This appears to happen 
when, due to a statistical fluctuation, the sample autocorrelation function p(t) has 
a sharp bend somewhere in its tail. 
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a a' 


r t' 


Energy 


T T 


T T 


E o = 


T T 


T 1 


E t = 2(J' + K) 


T 1 


T T 


E 2 = 2(J + 


T 1 


T 1 


£ 3 = 2(J + J') 



Table 1: Energies for a bond joining the spins (cr, r) and (cr', r'). These energies are 
invariant under the transformations cr, a 1 <->■ — cr, — cr' and r, r' <-> — r, — r'. On each 
bond we have added a constant energy J + J' + -K" to (2.1), in order to set E = 0. 



Point 


J = J' 


K 


2/ 


4-state Potts model 

ZF 

X2 

Ising model 


ilog3 ^ 0.274653 
ilogita, 0.302923 

4 A/2+V^-l 

-ilog (| - V2) « 0.344132 
ilog(l + V2) ^ 0.440687 


\\og?> ^ 0.274653 
ilog(l + y/2) « 0.220343 
1 b 6(5-3^2) _ 0.147920 

4 to 11-6V2 






1 

2 

« 0.735579 
1 



Table 2: Points of the self-dual curve of the symmetric AT model where our MC 
simulations were performed. The parameter y is defined in (2.9)/(2.11). We also 
include the values corresponding to the Ising model (DIs); the dynamic data corre- 
sponding to this point have been taken from Baillie and Coddington [7]. 



L 


X 


Ch (exact) 


T int,£ 


8 


41.392 ± 


0.008 


1.145559240 


2.589 ±0.005 


16 


139.58 


± 


0.04 


1.498704959 


3.258 ±0.005 


32 


470.12 


± 


0.20 


1.846767590 


4.016 ±0.005 


50 


1025.9 


± 


0.4 


2.069384825 


4.585 ±0.005 


64 


1581.4 


± 


0.5 


2.192211393 


4.899 ±0.010 


100 


3453.7 


± 


1.4 


2.413876309 


5.510±0.017 


128 


5319.2 


± 


2.4 


2.536331335 


5.874 ±0.016 


256 


17900 


± 


7.0 


2.879786255 


6.928 ±0.030 


512 


60185 


±28. 


3.222907954 


8.144 ±0.055 



Table 3: Data for the 2D Ising model. The susceptibility \ and the integrated 
autocorrelation time T int ,£ are taken from Ref. [7]. The value of the specific heat Ch 
is obtained from the exact formula of Ferdinand and Fisher [53]. 
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L 


MCS 


X 






Tmt,S 


T int,M 


2 


16 
32 
64 
128 
256 
512 
1024 


0.9 
1.9 
4.4 
2.9 
2.9 
2.9 
0.8 


141.41 ± 0.29 
474.23 ± 0.94 
1589. 67± 2.84 
5316.35 ± 16.49 
17771.34 ± 74.91 
59876.54 ± 334.67 
196872.43 ±3137.60 


5.027±0.027 
8.341 ±0.040 
13.937±0.060 
23.576±0.170 
39.876±0.388 
67.938 ±0.934 
120.270±4.183 


15.758 ± 0.056 
31.681 ± 0.100 
63.827± 0.172 
127.974 ± 0.575 
256.175 ± 1.523 
519.102 ± 4.078 
1020. 767±21. 680 


12.86 ± 0.24 
23.13 ± 0.40 
41.42 ± 0.62 
78.79 ± 2.01 
142.54 ± 4.90 
252.83 ±11.57 
534.44 ±67.68 


1 77 _|_ n 94 

22.75 ± 0.39 
40.34± 0.60 
76.19 ± 1.92 
136. 57± 4.59 
241.40 ±10.79 
505.27±62.22 



Table 4: Results of the MC simulations at the critical point of the 4-state Potts 
model. For each lattice size {L) we include the number of performed measurements 
(MCS) in units of 10 6 , the susceptibility (x), the specific heat {Ch), the second- 
moment correlation length (£), and the integrated autocorrelation times for the 
energy (ri ntj ,f) and the susceptibility (r int ^2). The quoted errors correspond to one 
standard deviation (i.e., confidence level 68%). 



L 


MCS 


Xlo 


XOT 




Car 


Cif,max 


Cif,min 


16 


0.9 


146.53 ± 0.24 


123.83 ± 0.25 


16.252 ±0.049 


13.751 ±0.039 


9.902 ±0.043 


0.4366 ±0.0014 


32 


0.9 


494. 14 ± 1.10 


395. 77± 1.08 


32.374±0.119 


27.289 ±0.094 


15.922 ±0.089 


0.4549 ±0.0025 


64 


0.9 


1668.42 ± 4.80 


1264.13 ± 4.47 


64.818 ±0.296 


54.410±0.229 


24.949±0.181 


0.4645 ±0.0044 


128 


0.9 


5665.88 ± 21.00 


4067.21 ± 18.60 


131.094±0.749 


109.699 ±0.577 


39.068±0.362 


0.4700 ±0.0050 


256 


0.9 


18925.33 ± 94.12 


12814.29 ± 78.80 


259.066 ±1.930 


216.909 ±1.488 


62.652 ±0.749 


0.4743 ±0.0086 


512 


1.9 


64118.12 ±274.65 


41094. 67±216. 67 


521.901 ±3.349 


436.890 ±2.548 


98.966 ±1.034 


0.4769 ±0.0085 



Table 5: Static data from the runs for the point ZF. For each lattice size L, we 
report the number of measurements (MCS) in units of 10 6 , the susceptibilities {\ w 
and x<jt), the second-moment correlation lengths {£ w and £ CTT ) and the maximum 
(C#,max) and the minimum {Ch, min) eigenvalues of the specific-heat matrix Ch- 



L 


MCS 






T int,Ml 


T int,Ml T 


16 


0.9 


9.43 ±0.15 


8.60 ±0.13 


9.23 ±0.15 


8.28 ±0.12 


32 


0.9 


16.00 ±0.33 


15.06 ±0.30 


15.39 ±0.31 


14.08 ±0.27 


64 


0.9 


26.40 ±0.70 


25.32 ±0.66 


24.98 ±0.65 


22.92 ±0.57 


128 


0.9 


44.97 ±1.56 


43.74 ±1.50 


41.25 ±1.37 


38.18 ±1.22 


256 


0.9 


76.35 ±3.45 


75.19 ±3.37 


70.28 ±3.05 


65.74 ±2.76 


512 


1.9 


119.02±4.62 


117.83±4.55 


110.03±4.11 


102.53 ±3.69 



Table 6: Autocorrelation times for the runs performed at the point ZF. For each 
lattice size L, we show the number of measurements (MCS) in units of 10 6 , the inte- 
grated autocorrelation times for the energies (7"^^ and Ti nt ,£^ T ) and the integrated 
autocorrelation times for the susceptibilities (T int ^2 and r int ^2 ). 
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L 


MCS 


Xlo 


XOT 




Car 


Cif,max 


Cif,min 


16 


0.9 


145.93 ± 0.18 


104.31 ± 


0.19 


15.815 ±0.036 


11.958±0.027 


8.793±0.030 


0.27326 ±0.00088 


32 


0.9 


493.30 ± 0.74 


319.01 ± 


0.72 


31.482 ±0.081 


23.779 ±0.059 


12.602±0.052 


0.28820 ±0.00090 


64 


0.9 


1660.56 ± 3.04 


970.81 ± 


2.63 


62.797±0.187 


47.488 ±0.134 


17.783 ±0.090 


0.29708 ±0.00140 


128 


0.9 


5595. 77± 12.41 


2958.76 ± 


9.63 


125.552 ±0.437 


94.974±0.312 


24.808±0.152 


0.30122±0. 00250 


256 


0.9 


18772.61 ± 48.29 


8951.11 ± 


33.94 


249.624±0.990 


188.436 ±0.696 


33.765 ±0.248 


0.30392 ±0.00316 


512 


0.9 


63352.87 ±189.64 


27341.35 ±120.27 


501.468±2.282 


378.482 ±1.609 


45.485 ±0.399 


0.30652 ±0.00462 



Table 7: The same as in Table 5 for the point X2. 



L 


MCS 








T int,Mlr 


16 


0.9 


6.405 ±0.085 


5.965 ±0.076 


6.172 ±0.081 


5.606 ±0.069 


32 


0.9 


9.326 ±0.148 


8.815 ±0.136 


8.706 ±0.134 


7.978 ±0.117 


64 


0.9 


13.686 ±0.264 


13.169 ±0.249 


12.372 ±0.227 


11.370 ±0.200 


128 


0.9 


20.338 ±0.476 


19.719 ±0.454 


17.798 ±0.389 


16.232 ±0.340 


256 


0.9 


27.775 ±0.758 


27.196 ±0.735 


23.777 ±0.601 


21.839 ±0.530 


512 


0.9 


39.559 ±1.288 


38.928 ±1.257 


32.474 ±0.957 


29.712 ±0.839 



Table 8: The same as in Table 6 for the point X2. 



Ratio 


4-state Potts model 
numerical exact 


ZF model 
numerical exact 


X2 model 
numerical exact 


Ising model 
numerical exact 


llv 

Hv 

a/v 


1.744 ±0.001 7/4 
1.744 ±0.001 7/4 
0.768 ± 0.009 lxlog- 3/2 
0.876 ±0.012 > 1 x log _3/2 


1.750 ± 0.004 7/4 
1.668 ± 0.005 5/3 
0.663 ± 0.006 2/3 
0.740 ±0.010 > 2/3 


1.751 ±0.001 7/4 
1.605 ± 0.001 1.6045 
0.438 ± 0.008 0.4183 
0.477 ± 0.028 > 0.4183 


1.7501 ± 0.0002 7/4 
1/2 
log 

0.240 ± 0.004 > log 



Table 9: Ratios of static critical exponents and dynamic critical exponent (zi nt ,£) 
coming from the power-law fits of the results contained in Tables 4-8. For the 
Ising model we include the fits to the dynamical data reported in Ref. [7]. For each 
model we present two columns, one with the MC results (the left one) and the other 
with the exact known results (the right one). The errors represent one standard 
deviation (i.e., confidence level of 68%). The symbol "1 X log -3 / 2 " means that the 
leading term of the specific heat for the 4-state Potts model behaves like L log -3 / 2 L. 
Likewise, the symbol "log" means that the leading term of the specific heat for the 
Ising model is logL. 
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L 


alv 
/ 


C H 


T 

~ L 1 


J 

x 2 








7"int,£ 


x 2 






16 


0.749 ± 0.003 


13.82 


(DF= 


5, 


level= 


2%) 


0.867 ±0.009 


4.32 


(DF= 


5, 


level= 


50%) 


32 


0.756 ± 0.004 


5.66 


(DF= 


4, 


level= 


23%) 


0.876 ± 0.012 


3.16 


(DF= 


= 4, 


level= 


53%) 


64 


0.762 ± 0.005 


1.84 


(DF= 


3, 


level= 


61%) 


0.887± 0.017 


2.26 


(DF= 


3, 


level= 


52%) 


128 


0.768 ± 0.009 


1.40 


(DF= 


= 2, 


level= 


50%) 


0.861 ±0.033 


1.31 


(DF= 


2, 


level= 


52%) 


256 


0.781 ±0.019 


0.71 


(DF= 


1, 


level= 


40%) 


0.880 ±0.067 


1.20 


(DF= 


1, 


level= 


27%) 


512 


0.824 ± 0.054 


0.00 


(DF= 


o, 


level= 


100%) 


1.080 ±0.194 


0.00 


(DF= 


o, 


level= 


100%) 



Table 10: Estimates of ajv and z- in t t £ for the 4-state Potts model at criticality 
as a function of the points involved in the fit (L > L m i n ). Errors represent one 
standard deviation, DF stands for the number of degrees of freedom and "level" is 
the confidence level of the fit (i.e., the probability that % 2 would equal or exceed 
the observed value, assuming that the underlying statistical model is correct). The 
preferred fits are marked with boldface. 







C H ~L 


a > v log 


p-3/2 


L 






C H 


~ L\og 


-p 


L 






ct/v 






X 


2 




P 








x 2 




16 


1.118 ± 0.003 


236.01 


(DF= 


5, 


level= 


5 x 10" 49 ) 


1.008 ±0.011 


86.06 


(DF= 


5, 


level= 


5 x 10" 17 ) 


32 


1.086 ±0.004 


54.23 


(DF= 


4, 


level= 


5 x 10" 11 ) 


1.102 ± 0.016 


24.02 


(DF= 


4, 


level= 


0.008%) 


64 


1.062 ±0.005 


9.47 


(DF= 


3, 


level= 


2%) 


1.185 ±0.025 


5.33 


(DF= 


3, 


level= 


15%) 


128 


1.039 ± 0.009 


0.45 


(DF= 


= 2, 


level= 


80%) 


1.286 ± 0.052 


0.34 


(DF= 


: 2 


, level= 


84%) 


256 


1.030 ±0.019 


0.19 


(DF= 


1, 


level= 


66%) 


1.320 ± 0.115 


0.23 


(DF= 


1, 


level= 


63%) 


512 


1.052 ±0.054 


0.00 


(DF= 


o, 


level= 


100%) 


1.158 ±0.355 


0.00 


(DF= 


o, 


level= 


100%) 



Table 11: Results of the weighted least-squares fits for the specific heat of the 4-state 
Potts model at criticality to a function Ch = AL a l v log -3 / 2 L (first column) and to 
a function Ch = A 1 L\og~ v L (second column). The exact results are ajv = 1 and 
p = | respectively [18,16,17]. For each fit we show the % 2 , the number of degrees of 
freedom (DF) and the confidence level ("level"). 



L 


ZF 


X2 


16 


-0.86034 ± 0.00065 


-0.77163 ± 0.00052 


32 


-0.88553 ± 0.00059 


-0.79294 ± 0.00042 


64 


-0.89995 ± 0.00052 


-0.80885 ± 0.00037 


128 


-0.90944 ± 0.00039 


-0.81835 ± 0.00034 


256 


-0.91529 ± 0.00035 


-0.82414 ± 0.00031 


512 


-0.91858 ± 0.00021 


-0.82828 ± 0.00028 


oo 


-0.92388 (exact) 


-0.83771 (exact) 



Table 12: Eigenvectors of the specific-heat matrix Ch- They are parametrized such 
that Wmi n = (1, a) corresponds to the smaller eigenvalue CH,min and w max = (a, — 1) 
to the largest eigenvalue CH,ma,x- For the models ZF and X2 we give the measured 
values of the parameter function of the lattice size L. The bottom row 

(L = oo) shows the theoretically predicted infinite- volume value of a taken from 
the (4.26) [see text]. 



50 







C H 


max ^ 


L a l 


V 






r int,£^ 










a/v 








x 2 




z int,£ w 






x 2 






16 


0.484 ± 0.002 


76.78 


(DF= 


4, 


level= 


8 x 10" 16 ) 


0.534 ±0.007 


5.16 


(DF= 


4, 


level= 


27%) 


32 


0.469 ± 0.003 


26.22 


(DF= 


3, 


level= 


9 x 10" 6 %) 


0.527 ± 0.011 


4.29 


(DF= 


3, 


level= 


23%) 


64 


0.455 ± 0.004 


7.50 


(DF= 


2, 


level= 


2%) 


0.509 ±0.016 


2.38 


(DF= 


2, 


level= 


30%) 


128 


0.438 ± 0.008 


0.32 


(DF= 


= 1, 


level= 


57%) 


0.477 ± 0.028 


0.39 


(DF= 


= 1, 


level= 


53%) 


256 


0.430 ± 0.017 


0.00 


(DF= 


o, 


level= 


100%) 


0.510 ±0.061 


0.00 


(DF= 


o, 


level= 


100%) 



Table 13: Estimates of a/v and ^mt,^ for the point X2 as a function of the points 
involved in the fit (L > L m i n ). We also show the % 2 , the number of degrees of 
freedom (DF) and the confidence level ("level") of each fit. 



Point 


Type 




x* (= x* = x* T ) 


x 2 


Potts 
Potts 


c 

L 


128 
16 


1.002 ±0.003 
1.023 ±0.007 


2.56 (DF= 3, level= 46%) 
1.94 (DF= 5, level= 86%) 




Point 


Type 




< 


x 2 


ZF 


c 


16 


1.015 ±0.002 


4.12 (DF= 5, level= 53%) 


X2 
X2 
X2 


c 

L 

* = i 


64 
16 
16 


0.980 ±0.002 
0.965 ±0.006 
0.975 ±0.003 


1.80 (DF= 3, level= 61%) 
1.15 (DF= 4, level= 89%) 
1.08 (DF= 4, level= 80%) 



Point 


Type 






x 2 


ZF 


c 


16 


0.852 ±0.002 


2.26 (DF= 4, level= 69%) 


X2 
X2 
X2 


c 

L 

A = | 


256 
16 
16 


0.737 ±0.002 
0.729 ±0.004 
0.736 ±0.002 


0.57 (DF= 1, level= 45%) 

2.08 (DF= 4, level= 72%) 

2.09 (DF= 4, level= 72%) 



Table 14: Estimates of x* for the three points of the AT self-dual curve considered in 
this paper. For each point we present the result of the least-square fit to a constant 
(C) or to a function of the type (4.44) (L). For the X2 model we also include the fit 
to constant plus corrections of order A = |. We also include the % 2 , the number of 
degrees of freedom (DF) and the confidence level ("level"). 
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X* 


x(L) 


= X* 


x 2 






x(L) = x* 

X* 


+ A/\og2L 

x 2 




16 


0.9943 ± 0.0015 


19.00 


(DF= 


6, 


level= 


0.4%) 


1.0228 ±0.0071 


1.95 


(DF= 


= 5, 


level= 


86%) 


32 


0.9965 ± 0.0017 


10.19 


(DF= 


5, 


level= 


7%) 


1.0274 ± 0.0107 


1.62 


(DF= 


4, 


level= 


80%) 


64 


0.9993 ± 0.0021 


3.98 


(DF= 


4, 


level= 


41%) 


1.0247 ± 0.0165 


1.58 


(DF= 


3, 


level= 


66%) 


128 


1.0023 ± 0.0032 


2.56 


(DF= 


= 3, 


level= 


46%) 


1.0380 ± 0.0326 


1.35 


(DF= 


2, 


level= 


51%) 


256 


1.0050 ± 0.0046 


1.92 


(DF= 


2, 


level= 


38%) 


1.0683 ± 0.0713 


1.12 


(DF= 


1, 


level= 


29%) 


512 


1.0118 ± 0.0075 


0.57 


(DF= 


1, 


level= 


45%) 


0.8436 ± 0.2235 


0.00 


(DF= 


o, 


level= 


100%) 


1024 


0.9968 ± 0.0211 


0.00 


(DF= 


o, 


level= 


100%) 

















x(L) = x 

X* 


* + a/Vl 
x 2 






X 

X* 


[L) = 


x* + A/L 

x 2 






16 


1.0098 ± 0.0041 


1.94 


(DF= 5, 


level= 


85%) 


1.0023 ± 0.0026 


3.38 


(DF= 5, 


level= 


64%) 


32 


1.0118 ± 0.0055 


1.65 


(DF= 4, 


level= 


80%) 


1.0051 ± 0.0035 


1.92 


(DF= 4, 


level= 


75%) 


64 


1.0109 ± 0.0078 


1.62 


(DF= 3, 


level= 


65%) 


1.0060 ± 0.0050 


1.86 


(DF= 3, 


level= 


60%) 


128 


1.0163 ± 0.0132 


1.36 


(DF= 2, 


level= 


51%) 


1.0100 ± 0.0080 


1.44 


(DF= 2, 


level= 


49%) 


256 


1.0283 ± 0.0257 


1.07 


(DF= 1, 


level= 


30%) 


1.0189 ± 0.0147 


0.93 


(DF= 1, 


level= 


33%) 


512 


0.9557 ± 0.0745 


0.00 


(DF= 0, 


level= 


100%) 


0.9798 ± 0.0431 


0.00 


(DF= 0, 


level= 


100%) 



Table 15: Estimates of x(L) as a function of the number of points included in the 
fit (L > L m i n ) for the 4-state Potts model at criticality. We show the fit of x(L) 
to a pure constant, as well as to a constant plus some decreasing functions. We 
also include the % 2 , the number of degrees of freedom (DF) and the confidence level 
("level") for each fit. 
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Ansatz 


Ising 


X2 


ZF 


4-state Potts 


ALP 




p = 0.060 ± 0.004 


p= 0.051 ±0.009 


p= 0.077 ±0.012 


p = 0.118 ± 0.012 






x 2 = 


1.02 (2 DF, 60%) 


x 2 = 


1.28 (4 DF, 86%) 


x 2 = 


1.23 (4 DF, 87%) 


x 2 = 


1.39 (5 DF, 93%) 








= 100 




= 16 




= 16 




= 16 


Alo 


y p L 


p = 0.315 ± 0.020 


p= 0.263 ±0.061 


p= 0.311 ±0.049 


p = 0.543 ±0.073 






x 2 = 


2.59 (2 DF, 27%) 


x 2 = 


0.63 (3 DF, 89%) 


x 2 = 


1.15 (4 DF, 89%) 


x 2 = 


1.92 (4 DF, 75%) 








= 100 




= 32 




= 16 




= 32 


A + 


BlogL 


\ 2 = 

A 


1.43 (2 DF, 49%) 


\ 2 = 

A 


1.43 (4 DF, 84%) 


\ 2 = 

A 


1.03 (4 DF, 91%) 


\ 2 = 

A 


1.98 (5 DF, 85%) 








= 100 




= 16 




= 16 




= 16 


A 




x 2 = 


36.74 (1 DF, 10" 9 ) 


x 2 = 


1.46 (2 DF, 48%) 


x 2 = 


0.77 (3 DF, 68%) 


x 2 = 


1.56 (2 DF, 46%) 








= 256 




= 128 




= 128 




= 256 


A + 


a 

logL 


x 2 = 


4.00 (1 DF, 5%) 


x 2 = 


0.90 (3 DF, 76%) 


x 2 = 


0.82 (3 DF, 85%) 


x 2 = 


0.60 (3 DF, 75%) 






= 128 




= 32 




= 32 




= 64 


A + 


a 


x 2 = 


2.57 (2 DF, 28%) 


x 2 = 


0.63 (3 DF, 89%) 


x 2 = 


1.01 (4 DF, 91%) 


x 2 = 


0.77 (3 DF, 86%) 






= 100 




= 32 




= 16 




= 64 


A + 


a 

L 1 / 4 


x 2 = 


4.04 (2 DF, 13%) 


x 2 = 


0.67 (3 DF, 88%) 


x 2 = 


1.33 (4 DF, 86%) 


x 2 = 


0.66 (3 DF, 88%) 






= 100 




= 32 




= 16 




= 64 


A + 


a 

VI 


x 2 = 


5.31 (1 DF, 2%) 


x 2 = 


0.99 (3 DF, 80%) 


x 2 = 


0.83 (3 DF, 84%) 


x 2 = 


0.64 (3 DF, 89%) 






= 128 




= 32 




= 32 




= 64 


A + 


a 

L 


x 2 = 


10.67 (1 DF, 0.1%) 


x 2 = 


0.64 (2 DF, 42%) 


x 2 = 


0.18 (2 DF, 91%) 


x 2 = 


1.11 (3 DF, 77%) 








= 128 




= 64 




= 64 




= 64 



Table 16: Results of fitting the ratio T- m t,e/CH for different Ansatze (power-law, 
logarithmic, and bounded) for all the models. We only show the "best" fits for each 
case. For each of them we give the value of the % 2 , the number of degrees of freedom 
(DF), the confidence level and the L m i n used. For the first two power-law fits we 
also give the estimate of that power. 





7"int,£ ~ 

P 


- L log" 


-PL 

x 2 






Tmt,£/C H ~ lc 

P 


gPL 

x 2 






16 


0.538 ± 0.038 


3.35 


(DF= 


5, 


level= 


65%) 


0.472 ±0.049 


3.64 


(DF= 


5, 


level= 


60%) 


32 


0.563 ± 0.056 


3.00 


(DF= 


4, 


level= 


56%) 


0.543 ±0.073 


1.92 


(DF= 


: 4, 


level= 


75%) 


64 


0.558 ± 0.087 


3.00 


(DF= 


3, 


level= 


39%) 


0.630 ±0.112 


0.85 


(DF= 


3, 


level= 


84%) 


128 


0.776 ±0.180 


1.09 


(DF= 


= 2, 


level= 


58%) 


0.515 ±0.232 


0.52 


(DF= 


2, 


level= 


77%) 


256 


0.736 ± 0.402 


1.08 


(DF= 


1, 


level= 


30%) 


0.593 ±0.518 


0.50 


(DF= 


1, 


level= 


48%) 


512 


-0.525 ± 1.280 


0.00 


(DF= 


o, 


level= 


100%) 


1.683 ± 1.632 


0.00 


(DF= 


o, 


level= 


100%) 



Table 17: Results for the 4-state Potts model at criticality of fitting the autocor- 
relation time T int ,£ to a function AL log~ p (first column) and of fitting the ratio 
Cff/f"mt,£ to a function Blog~ p L (second column). In both cases, we show the % 2 , 
the number of degrees of freedom (DF) and the confidence level ("level"). 
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train 


L = 16 

tmax = 52 


L = 32 

tmax = 92 


L = 64 

tmax = 164 


L = 128 

tmax = 320 


L = 256 

tmax = 560 


L = 512 

tmax = 1012 


1 


T exp>£ = 12.32 ± 0.07 
X 2 = 312.84, DF=50 
level=10- 39 


T exp>£ = 20.56 ±0.10 
X 2 = 1315.1, DF=90 
level=10- 216 


Texp,£ =35.21 ±0.14 
X 2 = 5411.0, DF=162 
level=6 X 10" 1020 


Tex P ,£ = 65.37 ±0.41 
X 2 = 7266.9, DF=318 
level=2 X 10" 1296 


r expj £ = 125.35 ± 0.95 
X 2 = 11584, DF=558 
level=2 X 10" 2030 




0-5r; nti£ 


Texp,£ = 13.43 ± 0.12 
X 2 = 48.84, DF=45 
level=32% 


Tex P ,£ = 23.62 ± 0.21 
X 2 = 111.17, DF=79 
level=l% 


Tex P ,£ = 42.99 ± 0.32 
X 2 = 148.32, DF=142 
level=34% 


NO 

CONVERGENCE 


NO 

CONVERGENCE 


NO 

CONVERGENCE 


Tint,£ 


Texp,£ = 13.69 ± 0.21 
X 2 = 31.24, DF=38 
level=77% 


Tex P ,£ = 24.33 ± 0.34 
X 2 = 82.06, DF=68 
level=12% 


Tex P ,£ =45.01 ±0.54 
X 2 = 130.94, DF=122 
level=27% 


Texp,£ = 100.1 ± 2.0 
X 2 = 193.13, DF=239 
level=99% 


NO 

CONVERGENCE 


Texp,£ = 334 ± 10 
X 2 = 96.08, DF=758 
level=100% 


l-5r; nti£ 


Tex P ,£ = 13.62 ± 0.32 
X 2 = 25.66, DF=32 
level=77% 


r exp> £ = 25.72 ±0.61 
X 2 = 57.99, DF=56 
level=40% 


Tex P ,£ = 45.99 ±0.91 
X 2 = 100.26, DF=101 
level=50% 


T-exp,£ = 105.6 ± 3.6 
X 2 = 113.74, DF=199 
level=100% 


Tex P ,£ = 185.3 ±8.3 
X 2 = 86.47, DF=349 
level=100% 


Texp,£ = 324 ±14 
X 2 = 52.06, DF=631 
level=100% 


2T; ntj£ 


Tex P ,£ = 13.57 ±0.55 
X 2 = 17.98, DF=25 
level=84% 


Texp,£ = 26.7 ±1.0 
X 2 = 40.33, DF=45 
level=67% 


Texp,£ = 51.0 ±1.5 
X 2 = 91.85, DF=122 
level=98% 


Texp,£ = 112.8 ±6.4 
X 2 = 65.00, DF=159 
level=100% 


r exp> £ = 232 ± 19 
X 2 = 43.11, DF=279 
level=100% 


Tex P ,£ = 329 ±24 
X 2 = 31.42, DF=505 
level=100% 



Table 18: Values of r exPj £ for the 4-state Potts model, using t max = 4Ti nt) £ (actually, 
the nearest integer to 4ri ntj ,f), and various values of t m i n . In each case we present the 
estimate of r exPj £, the % 2 value of the fit, the number of degrees of freedom (DF) and 
the confidence level. For each lattice size L we have used as i m ,- n the integer nearest 
to the value shown in the first column. Those fits where our self-consistent process 
did not converge are marked with "NO CONVERGENCE". The fit for L = 512 
and t m i n = 1 has not been performed, as it is very memory-consuming and it is 
expected to be rather poor. 



tmin 


L = 16 

tmax = 39 


L = 32 

tmax = 69 


L = 64 

tmax = 123 


L = 128 

tmax = 240 


L = 256 

tmax = 420 


L = 512 

tmax = 759 


i 


r exp> £ = 12.31 ±0.07 
X 2 = 301.70, DF=37 
level=3 X 10~ 43 


Texp,£ = 20.48 ± 0.10 
X 2 = 1239.0, DF=67 
level=3 X 10~ 215 


Tex P ,£ = 35.08 ±0.14 
X 2 = 5269.0, DF=121 
level=2 X 10" 1020 


Tex P ,£ = 63.24 ±0.40 
X 2 = 6390.1, DF=238 
level=2 X 10" 1169 


Tex P ,£ = 114.48 ±0.93 
X 2 = 8987.9, DF=418 
level=5 X 10" 1586 




0.5r; nti £ 


Texp,£ = 13.40 ±0.12 
X 2 = 42.82, DF=32 
level=10% 


Tex P ,£ = 23.44 ± 0.20 
X 2 = 75.97, DF=56 
level=4% 


Tex P ,£ = 42.70 ± 0.32 
X 2 = 107.2, DF=101 
level=32% 


Texp,£ =89.5 ±1.1 

X 2 = 404.50, DF=199 
level=4 X 10~ 16 


NO 

CONVERGENCE 


NO 

CONVERGENCE 


T"int,£ 


Texp,£ = 13.63 ± 0.21 
X 2 = 25.55, DF=25 
level=43% 


r exp> £ = 23.95 ±0.34 
X 2 = 50.84, DF=45 
level=25% 


Tex P ,£ = 44.28 ±0.54 
X 2 = 85.97, DF=81 
level=33% 


Texp,£ =94.0 ±2.0 
X 2 = 156.5, DF=159 
level=54% 


Tex P ,£ = 163.9 ±4.6 
X 2 = 125.8, DF=279 
level=100% 


Texp,£ = 329 ±12 
X 2 = 77.82, DF=503 
level=100% 


l-5r; nti £ 


Tex P ,£ = 13.50 ±0.32 
X 2 = 20.41, DF=19 
level=37% 


r exp> £ = 24.88 ±0.59 
X 2 = 31.04, DF=33 
level=56% 


Tex P ,£ = 44.46 ± 0.88 
X 2 = 58.50, DF=60 
level=53% 


Texp,£ = 97.0 ±3.5 
X 2 = 88.59, DF=119 
level=98% 


Texp,£ = 164.1 ± 7.6 
X 2 = 63.83, DF=209 
level=100% 


Texp,£ =314 ±18 

X 2 = 43.11, DF=376 
level=100% 


2T; ntj £ 


Tex P ,£ = 13.30 ±0.54 
X 2 = 12.18, DF=12 
level=43% 


r exp> £ = 25.10 ±0.94 
X 2 = 16.30, DF=22 
level=80% 


T-exp,£ = 45.4 ± 1.4 
X 2 = 43.52, DF=40 
level=32% 


Texp,£ = 101.8 ±6.0 
X 2 = 57.22, DF=79 
level=97% 


Texp,£ = 194 ± 17 

X 2 = 36.46, DF=139 
level=100% 


T-exp,£ = 313 ± 31 
X 2 = 27.20, DF=250 
level=100% 



Table 19: Values of r exPj £ for the 4-state Potts model, using t max = 3ri ntj £ (more 
precisely, the nearest integer to that value). Notation as in Table 18. 
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L = 16 

tmax = 26 


L = 32 

tmax = 37 


L = 64 

tmax = 55 


L = 128 

tmax = 81 


L = 256 

tmax = 111 


L = 512 

tmax — 158 


1 


Texp,£^ = 6.34 ± 0.03 
x 2 = 47.38, DF=24 
level= 0.03% 


T e xp,£„ =9.11 ±0.05 

X 2 = 150.26, DF=35 
level=3 X 10- 16 


Texp,£^ = 13.29 ±0.08 
X 2 = 265.08, DF=53 
level=6 X 10- 30 


T e xp,£„ = 18.94 ±0.13 
X 2 = 504.72, DF=79 
level=3 X 10- 63 


Texp,£^ = 26.45 ±0.21 
X 2 = 727.28, DF=109 
level=4 X 10~ 92 


r exp £ ^ = 36.86 ± 0.33 
X 2 = 958.10, DF=156 
level=2 X 10- 115 




Tex P ,£„ = 6.45 ± 0.04 
X 2 = 24.46, DF=22 
level=32% 


T-exp,£^ = 9.62 ± 0.08 
X 2 = 24.85, DF=31 
level=77% 


T-exp,£^ = 14.06 ±0.13 
X 2 = 80.84, DF=47 
level=0.2% 


T-exp,£^ = 21.03 ± 0.24 
X 2 = 95.18, DF=70 
level=2% 


T-exp,£^ = 30.61 ± 0.40 
X 2 = 152.19, DF=96 
level=2 X 10~ 4 


NO 

CONVERGENCE 


T int ,£ u 


t o — 6 58 + 07 
X 2 = 17.79, DF=19 
level=54% 


t o — 9 78 + 01? 
' exp,£ w — ajoiu.ii 

X 2 = 20.36, DF=27 

level=82% 


t o — 1445 + 93 
'exp,£ w — ±<±.<±o in u.zo 

X 2 = 69.09, DF=40 

level=0.3% 


t o — 91 74 + 40 
'exp,£ w — ii.141u.1u 

X 2 = 78.17, DF=60 

level=6% 


t o — 33 45 + 70 

6 X p o j^j — * * 

X 2 = 194.72, DF=72 
level=3 X 10~ 13 


t o — 47 8+1? 
'exp,£ w — <±*.oin ±.z 

X 2 = 104.12, DF=117 

level=90% 




Texp,£^ = 6.68 ± 0.13 
X 2 = 15.67, DF=15 
level=40% 


Texp,£^ = 9.96 ±0.21 
X 2 = 14.73, DF=22 
level=87% 


Texp,£^ = 15.23 + 0.41 
X 2 = 59.35, DF=33 
level=0.3% 


Tex P ,£^ = 23.14 ±0.73 
X 2 = 58.74, DF=49 
level=16% 


Texp,£^ = 34.2+1.2 
X 2 = 78.72, DF=68 
level=18% 


Texp,£^ = 48.0 ± 1.9 
X 2 = 74.83, DF=98 
level=96% 


2r int,£^ 


Texp,£^ = 6.64 ± 0.20 
X 2 = 9.41, DF=12 
level=67% 


Tex P ,£^ = 9.77 ±0.34 
X 2 = 16.30, DF=22 
level=80% 


Texp,£^ = 16.16 ±0.67 
X 2 = 44.94, DF=27 
level=2% 


Texp,£^ = 26.8 ±1.4 
X 2 = 37.44, DF=39 
level=54% 


Tex P ,£^ = 38.3 ± 2.3 
X 2 = 44.37, DF=54 
level=82% 


Texp,£^ = 48.1 ± 3.3 
X 2 = 43.41, DF=78 
level=100% 



Table 20: Values of T e ^s w for the X2 model, using t max = Ar-^t^^. Notation as in 
Table 18. 



train 


L = 16 

tmax = 19 


L = 32 

tmax = 28 


L = 64 

tmax = 41 


L = 128 

tmax = 61 


L = 256 

tmax = 83 


L = 512 

tmax = 119 


i 


Tex P ,£^ = 6.34 ± 0.03 
X 2 = 38.71, DF=17 
level= 0.2% 


Texp,£^ = 9.10 ±0.05 
X 2 = 141.89, DF=26 
level=6 X 10~ 18 


Tex P ,£^ = 13.23 ±0.08 
X 2 = 215.92, DF=39 
level=2 X 10- 26 


Texp,£^ = 18.72 ±0.13 
X 2 = 402.09, DF=59 
level=2 X 10- 52 


Tex P ,£^ = 26.09 ±0.21 
X 2 = 624.09, DF=81 
level=9 X 10- 85 


r expi£ ^ = 36.09 ± 0.33 
X 2 = 826.07, DF=117 
level=4 X 10- 107 


0.5r; nti£ ^ 


Tex P ,£^ = 6.45 ±0.04 
X 2 = 16.62, DF=15 
level=34% 


Texp,£^ = 9.61 ± 0.08 
X 2 = 18.96, DF=22 
level=65% 


Texp,£^ = 13.94 ±0.13 
X 2 = 44.77, DF=33 
level=8% 


Tex P ,£^ = 20.60 ±0.24 
X 2 = 40.85, DF=50 
level=82% 


Tex P ,£^ = 29.68 ±0.40 
X 2 = 98.71, DF=68 
level=0.9% 


Tex P ,£^ = 41.98 ±0.67 
X 2 = 86.30, DF=98 
level=79% 


T int,£^ 


Tex P ,£^ = 6.57 ±0.07 
X 2 = 10.77, DF=12 
level=55% 


Texp,£^ = 9.75 ±0.12 
X 2 = 14.96, DF=18 
level=66% 


Texp,£^ = 14.17 ±0.23 
X 2 = 37.19, DF=26 
level=7% 


Tex P ,£^ = 20.80 ±0.39 
X 2 = 28.85, DF=40 
level=90% 


Tex P ,£^ = 29.82 ±0.65 
X 2 = 79.36, DF=54 
level=l% 


Texp,£^ = 43.7 ± 1.2 
X 2 = 60.62, DF=78 
level=93% 




Texp,£^ = 6.64 ± 0.13 
X 2 = 8.90, DF=8 
level=35% 


Texp,£^ = 9.89 ±0.21 
X 2 = 9.82, DF=13 
level=71% 


Tex P ,£^ = 14.59 ±0.39 
X 2 = 28.89, DF=19 
level=7% 


Tex P ,£^ = 21.09 ±0.68 
X 2 = 19.59, DF=29 
level=91% 


r exp,£ w = 31.1 ± 1.1 
X 2 = 50.02, DF=40 
level=13% 


Texp,£^ = 42.3 ± 1.8 
X 2 = 46.00, DF=59 
level=89% 


2r int,£^ 


Tex P ,£^ = 6.56 ± 0.20 
X 2 = 2.58, DF=5 
level=76% 


Tex P ,£^ = 9.63 ±0.34 
X 2 = 8.18, DF=8 
level=42% 


Tex P ,£^ = 15.08 ±0.62 
X 2 = 20.69, DF=13 
level=8% 


Texp,£^ = 22.9 ±1.2 
X 2 = 10.32, DF=19 
level=94% 


Texp,£^ = 33.5 ±2.1 
X 2 = 27.71, DF=26 
level=37% 


Tex P ,£^ = 39.9 ± 2.8 
X 2 = 30.36, DF=39 
level=84% 



Table 21: Values of T e ^s w for the X2 model, using t max = Sr^t^^- Notation as in 
Table 18. 



tmin 


L = 16 

tmax = 38 


L = 32 

tmax = 64 


L = 64 

tmax = 106 


L = 128 

tmax = 180 


L = 256 

tmax = 306 


L = 512 

tmax = 476 


i 


Texp,£^ = 9.30 ± 0.05 
X 2 = 116.93, DF=36 
level=2 X 10~ 6 


Texp,£^ = 15.15 ±0.10 

X 2 = 358.79, DF=62 
level=2 X 10~ 23 


Tex P ,£^ = 24.39 ±0.19 
X 2 = 667.10, DF=104 
level=5 X 10~ 42 


Tex P ,£^ = 40.28 ±0.37 
X 2 = 1257.9, DF=178 
level=8 X 10~ 82 


Tex P ,£^ = 72.64 ±0.78 
X 2 = 2192.6, DF=304 
level=10- 282 


Texp,£^ = 113.5 ±1.1 

X 2 = 4274.1, DF=474 
level=5 X 10" 602 


0.5r; nti£ ^ 


Tex P ,£^ = 9.75 ± 0.08 
X 2 = 21.77, DF=32 
level=91% 


Texp,£^ = 16.64 ±0.17 
X 2 = 58.02, DF=55 
level=36% 


Tex P ,£^ = 28.08 ±0.36 
X 2 = 124.13, DF=92 
level=l% 


Tex P ,£^ = 52.25 ±0.83 
X 2 = 368.20, DF=156 
level=3 X 10- 19 


NO 

CONVERGENCE 


NO 

CONVERGENCE 


T int,£^ 


Texp,£^ = 9.83 ± 0.12 
X 2 = 21.06, DF=28 
level=82% 


Tex P ,£^ = 17.22 ±0.29 
X 2 = 48.94, DF=47 
level=40% 


Tex P ,£^ = 29.81 ± 0.63 
X 2 = 124.72, DF=79 
level=0.08% 


Texp,£^ = 57.9 ± 1.6 
X 2 = 137.91, DF=134 
level=39% 


Tex P ,£^ = 99.71 ± 3.2 
X 2 = 99.73, DF=229 
level=100% 


Texp,£^ = 147.1 ±4.2 
X 2 = 121.00, DF=356 
level=100% 




Texp,£^ = 9.85 ±0.21 
X 2 = 19.65, DF=23 
level=66% 


Tex P ,£^ = 17.80 ±0.49 
X 2 = 44.00, DF=39 
level=27% 


Texp,£^ = 31.2 ± 1.1 
X 2 = 85.57, DF=65 
level=4% 


Texp,£^ = 60.1 ± 2.6 
X 2 = 71.20, DF=111 
level=100% 


Texp,£^ = 104.7 ±5.6 
X 2 = 53.61, DF=190 
level=100% 


Texp,£^ = 160.5 ±7.1 
X 2 = 58.86, DF=296 
level=100% 


2r int,£^ 


Texp,£^ = 10.12 ±0.36 
X 2 = 17.26, DF=18 
level=51% 


Texp,£^ = 18.18 ±0.83 
X 2 = 39.46, DF=31 
level=14% 


Tex P ,£^ = 35.9 ± 2.3 
X 2 = 46.99, DF=52 
level=67% 


Tex P ,£^ =63.5 ±4.7 
X 2 = 39.02, DF=89 
level=100% 


Texp,£^ = 114 ±10 

X 2 =25.90, DF=152 
level=100% 


T-ex P ,£^ = 163 ± 13 
X 2 = 33.66, DF=237 
level=100% 



Table 22: Values of T e ^s w f° r the ZF model, using t max = Ar-^t^^. Notation as in 
Table 18. 
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L = 16 

^ Tftax — 28 


L = 32 

^ Tftax — 48 


L = 64 

~t ffidx — 79 


L = 128 

^ Tftax — 135 


L = 256 

^ Tftax — 229 


L = 512 

ifftax — 357 


1 


r expi£ ^ = 9.29 ± 0.05 
X 2 = 99.77, DF=26 
level=10- 10 


Texp,£^ = 15.10 ±0.10 

X 2 = 327.13, DF=46 
level=5 X 10~ 44 


r expj£ ^ = 24.14 ± 0.19 
X 2 = 579.33, DF=77 
level=5 X 10~ 78 


Tex P ,£^ = 38.93 ± 0.37 
X 2 = 1020.9, DF=133 
level=8 X 10~ 137 


Tex P ,£^ = 67.49 ± 0.79 
X 2 = 1692.7, DF=227 
level=3 X 10~ 222 


Texp,£^ = 107.5 ±1.1 

X 2 = 3642.5, DF=355 
level=10- 537 




Texp,£^ = 9.72 ± 0.08 
X 2 = 9.15, DF=22 
level=99% 


Texp,£^ = 16.55 ±0.17 
X 2 = 41.27, DF=39 
level=37% 


Tex P ,£^ = 27.54 ±0.36 
X 2 = 81.69, DF=65 
level=8% 


Tex P ,£^ = 50.57 ±0.86 
X 2 = 314.16, DF=111 
level=3 X 10~ 21 


NO 

CONVERGENCE 


Tex P ,£^ = 139.5 ±2.6 
X 2 = 213.48, DF=296 
level=100% 


T int,£^ 


Texp,£^ = 9.78 ± 0.12 
X 2 = 8.74, DF=18 
level=97% 


Texp,£^ = 17.01 ±0.29 
X 2 = 34.29, DF=31 
level=31% 


Tex P ,£^ = 28.17 ± 0.60 
X 2 = 69.80, DF=52 
level=5% 


r exp,£ w = 51.7 ± 1.5 

X 2 = 81.87, DF=89 
level=69% 


Tex P ,£^ = 93.9 ±3.5 
X 2 = 91.88, DF=152 
level=100% 


Tex P ,£^ = 137.6 ±4.2 
X 2 = 101.06, DF=237 
level=100% 




Texp,£^ = 9.72 ± 0.21 
X 2 = 7.37, DF=13 
level=88% 


Tex P ,£^ = 17.36 ± 0.48 
X 2 = 29.87, DF=23 
level=15% 


Texp,£^ = 28.8 ±1.1 

X 2 = 53.00, DF=38 
level=5% 


Texp,£^ = 51.8 ±2.5 
X 2 = 47.82, DF=66 
level=96% 


Tex P ,£^ = 93.9 ± 5.8 
X 2 = 43.53, DF=113 
level=100% 


r expi£ ^ = 148.0 ± 8.0 
X 2 = 50.63, DF=177 
level=100% 


2r int,£^ 


Tex P ,£^ = 9.84 ± 0.35 
X 2 = 5.54, DF=8 
level=70% 


Tex P ,£^ = 17.24 ±0.79 
X 2 = 26.92, DF=15 
level=3% 


Texp,£^ =31.6 ±2.0 
X 2 = 29.63, DF=25 
level=24% 


Tex P ,£^ = 52.3 ±4.3 
X 2 = 31.20, DF=44 
level=93% 


Texp,£^ = 100 ±10 
X 2 = 24.34, DF=75 
level=100% 


Texp,£^ = 145 ± 13 

X 2 = 32.14, DF=118 
level=100% 



Table 23: Values of T e ^s w for the ZF model, using t max = Sr^t^^- Notation as in 
Table 18. 



L 


4-state Potts model 

tmin tmax 7"int ,£ /^~exp } £ 




ZF model 

tmax 7"int ,£^ /^*exp,f w 


tmin 


X2 model 


16 


\ T mt,£ 


4r mt)f 


0.96 ±0.03 






0.97 ± 0.02 


2 T int,f^ 


4T m t,f^ 


0.99 ± 0.02 


32 


Tmt,S 


4r mt)f 


0.95 ±0.04 


2 T mt,£„ 


4-r mt) ^ 


0.96 ± 0.03 


2 T int,£^ 


4T m t,f^ 


0.97 ± 0.02 


64 


Tmt,S 


4r mt)f 


0.92 ±0.02 




3 Tint, e w 


0.96 ± 0.04 


2 T int,f^ 


3 T m t, £ w 


0.98 ± 0.03 


128 


Tmt,S 


3T m t,f 


0.84 ±0.04 


T mt,£„ 


3 Tint, e w 


0.87 ± 0.06 


^ T mt,£„ 


3 T m t, £ w 


0.99 ± 0.03 


256 


Tmt,S 


3T m t,f 


0.87 ±0.05 


T mt,£„ 


3 T m t, e w 


0.81 ±0.07 


^ T mt,£„ 


3 T m t, £ w 


0.94 ± 0.06 


512 


Tmt,S 


3T m t,f 


0.77 ±0.06 




3 T m t, e w 


0.85 ± 0.05 


2 T int,f^ 


3 T m t, £ w 


0.94 ± 0.05 



Table 24: Ratios T int ,£ / T eX p,£ for the three models considered in this paper. For each 
lattice size L we give the ratio and the interval [t m i n} t max ] used for its computation. 
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Figure 1: Phase diagram of the symmetric Ashkin-Teller model on the square lat- 
tice. The self-dual curve is B DIs P C. The solid curves represent second-order 
phase transitions, the dash-dotted ones the 4-state Potts-model subspace, and the 
dotted one the non-critical part of the self-dual curve. The Roman numerals desig- 
nate the different phases of the model (see text). 



57 



4 — State Potts Model Autocorrelation Function 

1 1^; i i i i I i i i i I i i i i I i i i i I 



H 0.1 
H 



- L=32 




L=16 J - 



0.01 







2 3 



t/r 



int.E 



4 



Figure 2: Autocorrelation function pssif) f° r the 4-state Potts model and L 
(□) and L = 32 (O), with the abscissa scaled by Ti ntj £ . The error bars are the s 
root of the diagonal terms of the covariance matrix (see Appendix B). 
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